diff options
Diffstat (limited to 'xicc/xlutfix.c')
| -rw-r--r-- | xicc/xlutfix.c | 1306 |
1 files changed, 1306 insertions, 0 deletions
diff --git a/xicc/xlutfix.c b/xicc/xlutfix.c new file mode 100644 index 0000000..6497688 --- /dev/null +++ b/xicc/xlutfix.c @@ -0,0 +1,1306 @@ +/* + * International Color Consortium color transform expanded support + * Set Lut table values and do auxiliary chanel interpolation continuity fixups. + * + * Author: Graeme W. Gill + * Date: 17/12/00 + * Version: 1.00 + * + * Copyright 2000 Graeme W. Gill + * All rights reserved. + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + * + */ + +/* + * This module provides additional xicc functionality + * for CMYK lut based profiles. + * + * This is essentially a test of one approach to fixing + * auxiliary parameter induced interpolation errors. + */ + +/* + * TTBD: + * + * Remove this code when the optimised separation code is working. + * + * Some of the error handling is crude. Shouldn't use + * error(), should return status. + * + */ + +/* Description: + + In all the clut based color systems, there are various + stages where the multi-dimenional profile functions are + resampled from one respresentation to another. As in all + sampling, aliasing may be an issue. The standard + methods for dealing with aliasing involve band limiting and + over-sampling. In dealing with color, anything other than + point sampling is often too slow to consider, meaning that + over sampling, or on-the-fly filtering is impractical. + + Band limiting the function to be sampled is therefore the + most practical approach, but there are still sever tradoffs. + For accurately representing the sampled characteristics + of a device, a high resolution grid, with band limited + sample points is desirable. 3 or 4 dimension grids however, + quickly consume memory, and generaly show an exponential + decline in access and manipulation speed with grid resolution. + To maintain accuracy therefore, the minimum grid resolution, + and the minimum level of filtering is often employed. + + The routines in this file are to deal with an aditional + subtlety when dealing with devices that have extra + degrees of freedom (ie. CMYK devices). In theory, the + aditional degrees of freedom can be set abitrarily, and + are often chosen to follow a "rule", designed to acheive + a goal such as minimising the amount of black used + in the highlights of bitonal devices (to minimise + "black dot" visibility), or to maximise black usage + for minimum ink costs, to resduce grey axis sensitivity + to the CMY values, to reduce the total ink loading, + or to pass through the inking values of a similar + input colorspace. In Argyll, these extra degrees of + freedom are refered to as auxiliary chanels. + + Because the final representation of the color correction + transform is often a multi-dimensional interpolation lookup + table (clut), there is an aditional hidden requirement for + any auxiliary input chanels, and that is that there be + a reasonable degree of interpolation continuity between + the sampled grid points. If this continuity requirement + is not met, then the accuracy of the interpolation within + each grid cell can be wildly inacurate, even though the + accuracy of the grid points themselves is good. + + For instance, if we have two grid points of a Lab->CMYK + interpolation grid: + + 1) 50 0 0 -> .0 .0 .0 .3 -> 50 0 0 + 2) 50 0 10 -> .2 .2 .4 .0 -> 50 0 10 + + Now if an input PCS value of 50 0 5 is used to + lookup the device values that should be used, a typical + interpolation might return: + + 50 0 5 -> .1 .1 .2 .15 -> 40 -5 6 + + This is a small change in PCS space, but bevcause the + two device points are at opposite extremes of the possible + auxliary locus for each point, the device values are + far appart in device space. The accuracy of the + device space interpolation is therefore not guaranteed + to be accurate, and might in this case, mean that + the device actually reproduces an unexpectedly inaccurate PCS value. + Even worse, at the gamut boundaries the locus shrinks to zero, + and particularly in the dark end of the gamut, there + may be a multitude of different Device values that overlap + at the gamut boundary, causing abrupt or even chaotic + device values at spacings well above the sampling spacing + of the interpolation grid being created. + + An additional challenge is that the locus of valid auxiliary + values may be discontinuous, (typically bifurcated), particularly + when an ink limit is imposed - the limit often removing a segment + of the auxiliary locus from the gamut. So ideally, a contiguous + auxiliary region needs to be mapped out, and any holes + patched over or removed from the gamut in a way that + doesn't introduce discontinuities. + + In Argyll, we want to maintain the freedom to set arbitrary + auxiliary rules, yet we need to avoid the gross loss of + accuracy abrupt transitions in auxiliary values can cause. + + The approach I've taken here involves a number of steps. The + first step sets up the clut in the usual maner, but records + various internal values for each point. In the second step, these + grid values are examined to locate cells which are "at risk" of + auxiliary interpolation errors. In the third step, the grid points + around the "at risk" cells have their auxiliary target values + adjusted to new target values, by using a simple smoothing filter + to reduce abrupt transitions. In the fourth step, new device values + are searched for, that have the same target PCS for the grid point, + but the smoothed auxiliary value. In cases where there is no scope + for meeting the new auxiliary target, because it is already at one + extreme of its possible locus for the target PCS, a tradoff is then + made between reduced target PCS accuracy, and an improved auxiliary + accuracy. In the final step, the new grid values replace the old + in the ICC. + +*/ + +#include "icc.h" +#include "numlib.h" +#include "xicc.h" + +/* NOTE:- that we only implement support for CMYK output here !!! */ + +#define CHECK_FUNCS /* Sanity check the callback functions */ +#define DO_STATS +#undef SAVE_TRACE /* Save the values returned by the clut callback function */ +#undef USE_TRACE /* Use the trace file instead of the clut callback function */ + +#define TRACENAME "D:/usr/argyll/xicc/xlutfix.xxx" /* So it will work in the debugger */ + +#define MAX_PASSES 7 +#define MAX_FILTERS 20 +#define THRESH 0.55 /* Fix Threshold, ratio of mean to maximum PCS point */ +#define MINTHRESH 2.0 /* Set minimum interp error threshold. Don't fix if below this */ +#define AUXWHT 3.0 /* Auxiliary tradeoff weight and increment */ +#undef WIDEFILTER /* Alter 4x4 neighborhood */ + +/* ========================================================== */ + +/* Return maximum difference */ +static double maxdiffn(int n, double *in1, double *in2) { + double tt, rv = 0.0; + int i; + for (i = 0; i < n; i++) { + if ((tt = fabs(in1[i] - in2[i])) > rv) + rv = tt; + } + return rv; +} + +/* Return absolute difference */ +static double absdiffn(int n, double *in1, double *in2) { + double tt, rv = 0.0; + int i; + for (i = 0; i < n; i++) { + tt = in1[i] - in2[i]; + rv += tt * tt; + } + return sqrt(rv); +} + +/* ========================================================== */ +/* Callback functions used by icc set_tables */ +/* ========================================================== */ + +/* Context for set_tables callbacks */ + +typedef struct { + int dofix; + void *cbctx; + void (*infunc)(void *cbctx, double *out, double *in); + void (*clutfunc)(void *cbctx, double *out, double *aux, double *auxr, double *pcs, double *in); + void (*clutpcsfunc)(void *cbctx, double *out, double *auxv, double *pcs); + void (*clutipcsfunc)(void *cbctx, double *pcs, double *olimit, double *auxv, double *in); + void (*outfunc)(void *cbctx, double *out, double *in); + + float *g; /* Base of grid */ + int res; /* Grid resolution */ + int fn; /* Number of floats in grid */ + int n; /* Number of entries in grid */ + int fesz; /* Entry size in floats */ + int fci[MXDI]; /* float increment for each input dimension into latice */ + int cmin[MXDI]; /* Fixup area bounding box minimum */ + int cmax[MXDI]; /* Fixup area bounding box maximum +1 */ + /* One float for flags */ + int din; /* Number of input (ie. grid) dimensions */ + int daux; /* Number of auxiliary dimensions */ + int dout; /* Number of output dimensions */ + int oauxr; /* Offset to start of aux range entries */ + int oauxv; /* Offset to start of aux value entries */ + int oauxvv; /* Offset to start of aux new value entries */ + int opcs; /* Offset to start of PCS value entries */ + int oout; /* Offset to start of output value entries */ + int nhi; /* Number of corners in an input grid cube */ + int *fhi; /* nhi grid cube corner offsets in floats */ + + /* Minimiser info */ + double m_auxw; /* Auxiliary error weighting factor (ie. 5 - 100) */ + double m_auxv[MAX_CHAN];/* Auxiliary target value */ + double m_pcs[3]; /* PCS target value */ + +#if defined(SAVE_TRACE) || defined(USE_TRACE) + FILE *tf; +#endif +} xifs; + +/* Macros to access flag values */ +#define XLF_FLAGV(fp) (*((unsigned int *)(fp))) +#define XLF_TOFIX 0x0001 /* Grid point to be fixed flag */ +#define XLF_UPDATE 0x0002 /* Grid point to be updated flag */ +#define XLF_HARDER 0x0004 /* Compromise PCS to improve result */ + +/* Functions to pass to icc settables() to setup icc Lut */ +/* Input table. input -> input' space. */ +static void xif_set_input(void *cntx, double *out, double *in) { + xifs *p = (xifs *)cntx; + + p->infunc(p->cbctx, out, in); +} + + +/* clut, input' -> output' space */ +static void xif_set_clut(void *cntx, double *out, double *in) { + xifs *p = (xifs *)cntx; + + if (p->dofix == 0) { /* No fixups */ + p->clutfunc(p->cbctx, out, NULL, NULL, NULL, in); + + } else if (p->dofix == 1) { /* First pass */ + int e, f; + float *fp, *ep; + double pcs[MAX_CHAN], auxv[MAX_CHAN], auxr[MAX_CHAN * 2]; + + /* the icclib set_tables() supplies us the grid indexes */ + /* as integer in the double locations at in[-e-1] */ + +#if defined(USE_TRACE) + if (fread(pcs, sizeof(double), 3, p->tf) != 3 + || fread(auxr, sizeof(double), 2 * p->daux, p->tf) != (2 * p->daux) + || fread(auxv, sizeof(double), p->daux, p->tf) != p->daux + || fread(out, sizeof(double), p->dout, p->tf) != p->dout) { + fprintf(stderr,"mark_cells: read of trace failed\n"); + exit(-1); + } +#else /* !USE_TRACE */ + p->clutfunc(p->cbctx, out, auxv, auxr, pcs, in); + +#if defined(SAVE_TRACE) + if (fwrite(pcs, sizeof(double), 3, p->tf) != 3 + || fwrite(auxr, sizeof(double), 2 * p->daux, p->tf) != (2 * p->daux) + || fwrite(auxv, sizeof(double), p->daux, p->tf) != p->daux + || fwrite(out, sizeof(double), p->dout, p->tf) != p->dout) { + fprintf(stderr,"mark_cells: write of trace failed\n"); + exit(-1); + } +#endif /* SAVE_TRACE */ +#endif /* !USE_TRACE */ + + /* Figure grid pointer to grid entry */ + for (fp = p->g, e = 0; e < p->din; e++) + fp += *((int *)&in[-e-1]) * p->fci[e]; + + XLF_FLAGV(fp) = 0; /* Clear flags */ + + ep = fp + p->opcs; + for (f = 0; f < 3; f++) /* Save PCS values */ + ep[f] = (float)pcs[f]; + + ep = fp + p->oauxr; + for (f = 0; f < (2 * p->daux); f++) /* Save auxiliary range values */ + ep[f] = (float)auxr[f]; + + ep = fp + p->oauxv; + for (f = 0; f < p->daux; f++) /* Save auxiliary values */ + ep[f] = (float)auxv[f]; + + ep = fp + p->oout; + for (f = 0; f < p->dout; f++) /* Save the output values */ + ep[f] = (float)out[f]; + + } else { /* Second pass */ + int e, f; + float *fp, *ep; + + /* Figure grid pointer to grid entry */ + for (fp = p->g, e = 0; e < p->din; e++) + fp += *((int *)&in[-e-1]) * p->fci[e]; + + ep = fp + p->oout; + for (f = 0; f < p->dout; f++) /* Return the fixed output values */ + out[f] = (double)ep[f]; + } +} + +/* output output' -> output space */ +static void xif_set_output(void *cntx, double *out, double *in) { + xifs *p = (xifs *)cntx; + + p->outfunc(p->cbctx, out, in); +} + +static int mark_cells(xifs *p); +static int filter_grid(xifs *p, int tharder); +static void fix_grid(xifs *p, double auxw); +static int comp_pcs(xifs *p, double auxw, double *auxrv, double *auxv, double *pcs, double *dev); + +/* Helper function to setup the three tables contents, and the underlying icc. */ +/* Only useful if there are auxiliary device output chanels to be set, */ +/* as this takes care of auxiliary interpolation continuity fixups. */ +int icxLut_set_tables_auxfix( +icmLut *p, /* Pointer to icmLut object */ +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) */ +double *inmin, double *inmax, /* Maximum range of inspace' values */ + /* (NULL = default) */ +void (*clutfunc)(void *cbctx, double *out, double *aux, double *auxr, double *pcs, double *in), + /* inspace' -> outspace' transfer function, also */ + /* return the target PCS and the (packed) auxiliary locus range, */ + /* as [min0, max0, min1, max1...], the auxiliary chosen. */ +void (*clutpcsfunc)(void *cbctx, double *out, double *aux, double *pcs), + /* PCS + aux_target -> outspace' transfer function */ +void (*clutipcsfunc)(void *cbctx, double *pcs, double *olimit, double *auxv, double *in), + /* outspace' -> PCS + auxv check 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) */ +) { + int rv, g, e, jj, kk; + double auxw; /* Auxiliary weight factor */ + xifs xcs; /* Our context structure */ + + /* Simply pass this on to the icc set_table() */ + xcs.dofix = 0; /* Assume we won't attempt fix */ + xcs.cbctx = cbctx; + xcs.infunc = infunc; + xcs.clutfunc = clutfunc; + xcs.clutpcsfunc = clutpcsfunc; + xcs.clutipcsfunc = clutipcsfunc; + xcs.outfunc = outfunc; + + if (outsig != icSigCmykData) { /* Don'y know how/if to fix this */ + rv = p->set_tables(p, + ICM_CLUT_SET_APXLS, + (void *)&xcs, + insig, outsig, + xif_set_input, + inmin, inmax, + xif_set_clut, + clutmin, clutmax, + xif_set_output); + + return rv; + } + +#ifdef CHECK_FUNCS + if (insig == icSigLabData) { + double in[3], out[MAX_CHAN]; + double aux[1], auxr[2], pcs[3]; + double out_check[MAX_CHAN]; + double apcs[3], pcs_check[3]; + + /* Pick a sample input value */ + in[0] = 50.0; in[1] = 0.0; in[2] = 0.0; + + /* Test the in->out function */ + clutfunc(cbctx, out, aux, auxr, pcs, in); + +printf("~1 %f %f %f -> pcs %f %f %f,\n auxr %f - %f, auxv %f, dev %f %f %f %f\n", + in[0], in[1], in[2], pcs[0], pcs[1], pcs[2], auxr[0], auxr[1], aux[0], + out[0], out[1], out[2], out[3]); + + /* Check that we get the same result for the pcs function */ + clutpcsfunc(cbctx, out_check, aux, pcs); + + if (maxdiffn(p->outputChan, out, out_check) > 1e-6) { + fprintf(stderr,"set_tables_auxfix: pcsfunc check failed\n"); +printf("~1 is %f %f %f %f, should be %f %f %f %f\n", +out_check[0], out_check[1], out_check[2], out_check[3], +out[0], out[1], out[2], out[3]); + } +printf("~1 PCS version gives %f %f %f %f\n", +out_check[0], out_check[1], out_check[2], out_check[3]); + + /* Checkout the reverse function */ + clutipcsfunc(cbctx, apcs, NULL, NULL, out); /* Device -> clipped PCS */ + +printf("~1 clipped PCS = %f %f %f\n", apcs[0], apcs[1], apcs[2]); + clutpcsfunc(cbctx, out_check, aux, apcs); /* clipped PCS -> Device */ + clutipcsfunc(cbctx, pcs_check, NULL, NULL, out_check); /* Device -> PCS */ +printf("~1 check PCS = %f %f %f\n", pcs_check[0], pcs_check[1], pcs_check[2]); + + if (maxdiffn(3, apcs, pcs_check) > 1e-5) { + fprintf(stderr,"set_tables_auxfix: ipcsfunc check failed\n"); + printf("~1 is %f %f %f, should be %f %f %f\n", + pcs_check[0], pcs_check[1], pcs_check[2], + pcs[0], pcs[1], pcs[2]); + } + + } else { + fprintf(stderr,"Sanity check of %s not implemented!\n", + icm2str(icmColorSpaceSignature,insig)); + } +#endif /* CHECK_FUNCS */ + + /* Allocate an array to hold all the results */ + xcs.res = p->clutPoints; + xcs.din = p->inputChan; + xcs.dout = p->outputChan; + xcs.daux = xcs.dout - 3; /* Number of auxiliary values */ + xcs.fesz = 1 + 3 + xcs.dout + 4 * xcs.daux; /* Entry size in floats */ + + /* Compute total number of entries, and offsets in each dimension */ + xcs.n = xcs.res; + xcs.fci[0] = xcs.fesz; + for (e = 1; e < xcs.din; e++) { + xcs.n *= xcs.res; + xcs.fci[e] = xcs.fci[e-1] * xcs.res; + } + xcs.fn = xcs.n * xcs.fesz; + +printf("~1 fci = %d %d %d\n", +xcs.fci[0], xcs.fci[1], xcs.fci[2]); + + /* Setup offset list to grid cube corners */ + xcs.nhi = 1 << xcs.din; + if ((xcs.fhi = (int *)malloc(sizeof(int) * xcs.nhi)) == NULL) { + sprintf(p->icp->err,"icxLut_set_tables: malloc() failed"); + return p->icp->errc = 2; + } + for (g = 0; g < xcs.nhi; g++) { + xcs.fhi[g] = 0; + for (e = 0; e < xcs.din; e++) { + if (g & (1 << e)) + xcs.fhi[g] += xcs.fci[e]; + } + + } +printf("~1 nhi = %dd\n",xcs.nhi); + + /* Offsets into each entry */ + xcs.opcs = 1; /* Allow 1 flag float */ + xcs.oout = xcs.opcs + 3; /* dpcs floats */ + xcs.oauxr = xcs.oout + xcs.dout; /* dout floats */ + xcs.oauxv = xcs.oauxr + xcs.daux * 2; /* 2 daux floats */ + xcs.oauxvv = xcs.oauxv + xcs.daux; /* daux floats */ + /* daux floats */ + +printf("~1 res %d, entry size = %d floats, total floats needed = %d\n",xcs.res,xcs.fesz,xcs.fn); + +printf("~1 opcs = %d, oout = %d, oauxr = %d, oauxv = %d\n", + xcs.opcs, xcs.oout, xcs.oauxr, xcs.oauxv); + + /* Allocate the grid */ + if ((xcs.g = (float *)malloc(sizeof(float) * xcs.fn)) == NULL) { + sprintf(p->icp->err,"icxLut_set_tables: malloc() failed"); + return p->icp->errc = 2; + } + +#if defined(SAVE_TRACE) || defined(USE_TRACE) + { + char *tname = TRACENAME; + +#if defined(SAVE_TRACE) + if ((xcs.tf = fopen(tname,"w")) == NULL) { +#else + if ((xcs.tf = fopen(tname,"r")) == NULL) { +#endif + fprintf(stderr,"mark_cells: Can't open file '%s'\n",tname); + exit(-1); + } +#if defined(O_BINARY) +#if defined(SAVE_TRACE) + xcs.tf = freopen(tname,"wb",xcs.tf); +#else + xcs.tf = freopen(tname,"rb",xcs.tf); +#endif +#endif + } +#endif /* SAVE_TRACE || USE_TRACE */ + +#ifdef NEVER +// ~9 check function +{ + int rv; + double auxv[1], rauxv[1]; + double pcs[3], rpcs[3]; + double dev[4]; + + auxv[0] = 0.5; + pcs[0] = 60.0; + pcs[1] = 0.0; + pcs[2] = 0.0; + + dev[0] = 0.5; + dev[1] = 0.1; + dev[2] = 0.1; + dev[3] = 0.1; + + rv = comp_pcs(&xcs, 20.0, NULL, auxv, pcs, dev); + + printf("~9 comp_pcs returned %d, device %f %f %f %f\n",rv, dev[0], dev[1], dev[2], dev[3]); + + xcs.clutipcsfunc(xcs.cbctx, rpcs, NULL, rauxv, dev); + + printf("~9 comp_pcs pcs %f %f %f, aux %f\n", rpcs[0], rpcs[1], rpcs[2], rauxv[0]); + + return 0; +} +#endif + +printf("~1 doing the first pass\n"); + /* First pass */ + xcs.dofix = 1; + rv = p->set_tables(p, + ICM_CLUT_SET_APXLS, + (void *)&xcs, + insig, outsig, + xif_set_input, + inmin, inmax, + xif_set_clut, + clutmin, clutmax, + xif_set_output); + +#if defined(SAVE_TRACE) || defined(USE_TRACE) + fclose(xcs.tf); +#endif /* SAVE_TRACE || USE_TRACE */ + + if (rv != 0) { + free(xcs.fhi); + free(xcs.g); + return rv; + } + +printf("~1 doing the fixups\n"); + + /* Try three passes */ + for(jj = 0, auxw = AUXWHT; jj < MAX_PASSES; jj++, auxw += AUXWHT) { + int lrv; + + /* Figure out which cells need fixing */ + rv = mark_cells(&xcs); +printf("~1 cells that need fixing = %d\n", rv); + + if (rv == 0) + break; + + /* Filter the grid values that need fixing */ +printf("~1 about to filter grid points\n"); + for (kk = 0, lrv = 0, rv = 1; kk < MAX_FILTERS && rv > 0 && rv != lrv; kk++) { + lrv = rv; + rv = filter_grid(&xcs, 1); + } + + if (rv == 0) + break; + + /* Lookup device values for grid points with changed auxiliary targets */ +printf("~1 about to fix grid points\n"); + fix_grid(&xcs, auxw); + }; + +rv = mark_cells(&xcs); +printf("~1 faulty cells remaining = %d\n", rv); + +printf("~1 updatding the icc\n"); + /* Second pass */ + xcs.dofix = 2; + rv = p->set_tables(p, + ICM_CLUT_SET_APXLS, + (void *)&xcs, + insig, outsig, + xif_set_input, + inmin, inmax, + xif_set_clut, + clutmin, clutmax, + xif_set_output); + + free(xcs.fhi); + free(xcs.g); + +printf("~1 done\n"); + + return rv; +} + + + + +/* ----------------------------------------- */ +/* Mark cells that need fixing */ +/* Return number of cells that need fixing */ +static int mark_cells(xifs *p) { + int e, f; + int coa[MAX_CHAN], ce; /* grid counter */ + int tcount = 0; +#ifdef DO_STATS + double aerr = 0.0; + double merr = 0.0; + double ccount = 0.0; +#endif + + /* Get ready to track fixup area bounding box */ + for (e = 0; e < p->din; e++) { + p->cmin[e] = 99999; + p->cmax[e] = -1; + } + + /* Init the grid counter */ + for (ce = 0; ce < p->din; ce++) + coa[ce] = 0; + ce = 0; + + /* Itterate throught the PCS clut grid cells */ + while (ce < p->din) { + int j, m; + float *gp; /* Grid pointer */ + float *ep, *fp; /* Temporary grid pointers */ + double wpcsd; /* Worst case PCS distance */ + double apcs[3]; /* Average PCS value */ + double aout[MAX_CHAN]; /* Average output value */ + double check[3]; /* Check PCS */ + double ier; /* Interpolation error */ + int markcell = 0; /* Mark the cell */ + + /* Compute base of cell pointer */ + gp = p->g; /* Grid pointer */ + for (e = 0; e < p->din; e++) + gp += coa[e] * p->fci[e]; + + /* - - - - - - - - - - - - - - - - - */ + /* Full surrounding Cell calculation */ + + /* Init averaging accumulators */ + for (j = 0; j < 3; j++) + apcs[j] = 0.0; + for (f = 0; f < p->dout; f++) + aout[f] = 0.0; + + /* For each corner of the PCS grid based at the current point, */ + /* average the PCS and Device values */ + for (m = 0; m < p->nhi; m++) { + double pcs[3]; + double dev[MAX_CHAN]; + + fp = gp + p->fhi[m]; + +//ep = fp + p->opcs; +//printf("Input PCS %f %f %f\n", ep[0], ep[1], ep[2]); + + ep = fp + p->oout; /* base of device values */ + for (f = 0; f < p->dout; f++) { + double v = (double)ep[f]; + dev[f] = v; + aout[f] += v; + } + + /* Device to clipped PCS */ + p->clutipcsfunc(p->cbctx, pcs, NULL, NULL, dev); + + for (j = 0; j < 3; j++) + apcs[j] += pcs[j]; + +//printf("Corner PCS %f %f %f -> ", pcs[0], pcs[1], pcs[2]); +//printf("%f %f %f %f\n", dev[0], dev[1], dev[2], dev[3]); + } + + for (j = 0; j < 3; j++) + apcs[j] /= (double)p->nhi; + + for (f = 0; f < p->dout; f++) + aout[f] /= (double)p->nhi; + + /* Compute worst case distance of PCS corners to average PCS */ + wpcsd = 0.0; + for (m = 0; m < p->nhi; m++) { + double ss; + + fp = gp + p->fhi[m] + p->opcs; + for (ss = 0.0, j = 0; j < 3; j++) { + double tt = (double)fp[j] - apcs[j]; + ss += tt * tt; + } + ss = sqrt(ss); + if (ss > wpcsd) + wpcsd = ss; + } + wpcsd *= THRESH; /* Set threshold as proportion of most distant corner */ + if (wpcsd < MINTHRESH) /* Set a minimum threshold */ + wpcsd = MINTHRESH; + +//printf("Average PCS %f %f %f, Average Device %f %f %f %f\n", +//apcs[0], apcs[1], apcs[2], aout[0], aout[1], aout[2], aout[3]); + + /* Average Device to PCS */ + p->clutipcsfunc(p->cbctx, check, NULL, NULL, aout); + +//printf("Check PCS %f %f %f\n", +//check[0], check[1], check[2]); + + /* Compute error in PCS vs. Device interpolation */ + ier = absdiffn(3, apcs, check); + +//printf("Average PCS %f %f %f, Check PCS %f %f %f, error %f\n", +//apcs[0], apcs[1], apcs[2], check[0], check[1], check[2], ier); + +#ifdef DO_STATS + aerr += ier; + ccount++; + if (ier > merr) + merr = ier; +#endif /* STATS */ + + if (ier > wpcsd) { + markcell = 1; +printf("~1 ier = %f, wpcsd = %f, Dev = %f %f %f %f\n", ier, wpcsd, aout[0], aout[1], aout[2], aout[3]); + } + + /* - - - - - - - - - - - - - */ + /* Point pair calculations */ + + if (markcell == 0) { /* Don't bother testing if already marked */ + + /* Get the base point values */ + ep = gp + p->oout; /* base of device values (assumes fhi[0] == 0) */ + + for (f = 0; f < p->dout; f++) + aout[f] = (double)ep[f]; + + p->clutipcsfunc(p->cbctx, apcs, NULL, NULL, aout); + + /* For each other corner of the PCS grid based at the */ + /* current point, compute the interpolation error */ + for (m = 1; m < p->nhi; m++) { + double pcs[3]; /* Average PCS */ + double dpcs[3]; /* PCS of averaged device */ + double dev[MAX_CHAN]; + + fp = gp + p->fhi[m]; /* Other point point */ + ep = fp + p->oout; /* base of device values */ + + for (f = 0; f < p->dout; f++) + dev[f] = (double)ep[f]; + + /* Device to clipped PCS */ + p->clutipcsfunc(p->cbctx, pcs, NULL, NULL, dev); + + for (j = 0; j < 3; j++) + pcs[j] = 0.5 * (pcs[j] + apcs[j]); /* PCS averaged value */ + + for (f = 0; f < p->dout; f++) + dev[f] = 0.5 * (aout[f] + dev[f]); /* Average device */ + + /* Average Device to PCS */ + p->clutipcsfunc(p->cbctx, dpcs, NULL, NULL, dev); + + wpcsd = THRESH * absdiffn(3, pcs, apcs); /* Threshold value */ + if (wpcsd < MINTHRESH) /* Set a minimum threshold */ + wpcsd = MINTHRESH; + + /* Compute error in PCS vs. Device interpolation */ + ier = absdiffn(3, pcs, dpcs); + +#ifdef DO_STATS + aerr += ier; + ccount++; + if (ier > merr) + merr = ier; +#endif /* STATS */ + if (ier > wpcsd) { /* Over threshold */ + markcell = 1; +printf("~1 ier = %f, wpcsd = %f, Dev = %f %f %f %f\n", ier, wpcsd, aout[0], aout[1], aout[2], aout[3]); + } + } + } + + if (markcell) { + +#ifndef WIDEFILTER + /* Mark the whole cube around this base point */ + tcount++; + /* Grid points that make up cell */ + for (m = 0; m < p->nhi; m++) { + float *fp = gp + p->fhi[m]; + XLF_FLAGV(fp) |= XLF_TOFIX; + } + for (e = 0; e < p->din; e++) { + if (coa[e] < p->cmin[e]) + p->cmin[e] = coa[e]; + if ((coa[e]+2) > p->cmax[e]) + p->cmax[e] = coa[e] + 2; + } +#else + int fo[MAX_CHAN], fe; /* region counter */ + + /* Mark the whole cube around this base point */ + tcount++; + /* Grid points one row beyond cell */ + for (fe = 0; fe < p->din; fe++) + fo[fe] = -1; /* Init the neighborhood counter */ + fe = 0; + + /* For each corner of the filter region */ + while (fe < p->din) { + float *fp = gp; + + for (e = 0; e < p->din; e++) { + int oo = fo[e] + coa[e]; + if (oo < 0 || oo >= p->res) + break; + if (oo < p->cmin[e]) + p->cmin[e] = oo; + if ((oo+1) > p->cmax[e]) + p->cmax[e] = oo + 1; + fp += fo[e] * p->fci[e]; + } + + if (e >= p->din) { /* Within the grid */ + XLF_FLAGV(fp) |= XLF_TOFIX; + } + + /* Increment the counter */ + for (fe = 0; fe < p->din; fe++) { + if (++fo[fe] < 3) + break; /* No carry */ + fo[fe] = -1; + } + } +#endif /* WIDEFILTER */ + } + + /* - - - - - - - - - - - - - - */ + + + /* Increment the main grid counter */ + for (ce = 0; ce < p->din; ce++) { + if (++coa[ce] < (p->res-1)) + break; /* No carry */ + coa[ce] = 0; + } + } + +#ifdef DO_STATS + aerr /= ccount; + + printf("~1Average interpolation error %f, maximum %f\n",aerr, merr); + printf("~1Number outside corner radius = %f%%\n",(double)tcount * 100.0/ccount); + printf("~1Bounding box is %d - %d, %d - %d, %d - %d\n", + p->cmin[0], p->cmax[0], p->cmin[1], p->cmax[1], p->cmin[2], p->cmax[2]); +#endif /* STATS */ + + return tcount; +} + + +/* ----------------------------------------- */ +/* Do one filter pass of grid aux values that need fixing */ +/* If tharder is set, don't clamp new aux targets, but signal trading PCS for aux. */ +/* Return the number of grid points that have a new aux target */ +static int filter_grid(xifs *p, int tharder) { + int e, f; + int coa[MAX_CHAN], ce; /* grid counter */ + int tcount = 0; + + /* Init the grid counter */ + for (ce = 0; ce < p->din; ce++) + coa[ce] = p->cmin[ce]; + ce = 0; + + /* Itterate throught the PCS clut grid cells, */ + /* computing new auxiliary values */ + while (ce < p->din) { + float *gp, *fp; /* Grid pointer */ + int fo[MAX_CHAN], fe; /* filter counter */ + double faux[MAX_CHAN]; /* Filtered auxiliary value */ + double nfv; /* Number of values in filter value */ + + /* Compute base of cell pointer */ + gp = p->g; /* Grid pointer */ + for (e = 0; e < p->din; e++) + gp += coa[e] * p->fci[e]; + + if ((XLF_FLAGV(gp) & XLF_TOFIX) != 0) { + + for (f = 0; f < p->daux; f++) + faux[f] = 0.0; + nfv = 0.0; + + /* Init the neighborhood counter */ + for (fe = 0; fe < p->din; fe++) + fo[fe] = -1; + fe = 0; + + /* For each corner of the filter region, */ + /* compute a filter kernel value */ + while (fe < p->din) { + + fp = gp + p->oauxv; + for (e = 0; e < p->din; e++) { + int oo = coa[e] + fo[e]; + if (oo < 0 || oo >= p->res) + break; + fp += fo[e] * p->fci[e]; + } + + if (e >= p->din) { /* Add this neighborhood value */ + for (f = 0; f < p->daux; f++) + faux[f] += fp[f]; + nfv++; + } + + /* Increment the counter */ + for (fe = 0; fe < p->din; fe++) { + if (++fo[fe] < 2) + break; /* No carry */ + fo[fe] = -1; + } + } + + /* Compute average value, and save it */ + fp = gp + p->oauxvv; + for (f = 0; f < p->daux; f++) + fp[f] = (float)(faux[f] / nfv); + + } + + /* Increment the counter */ + for (ce = 0; ce < p->din; ce++) { + if (++coa[ce] < p->cmax[ce]) + break; /* No carry */ + coa[ce] = p->cmin[ce]; + } + } + + /* Clip the new values to the valid range, and put them into place */ + + /* Init the grid counter */ + for (ce = 0; ce < p->din; ce++) + coa[ce] = p->cmin[ce]; + ce = 0; + + /* Itterate throught the PCS clut grid cells, */ + /* computing new auxiliary values */ + while (ce < p->din) { + float *gp, *dp, *sp, *rp; /* Grid pointer */ + + /* Compute base of cell pointer */ + gp = p->g; /* Grid pointer */ + for (e = 0; e < p->din; e++) + gp += coa[e] * p->fci[e]; + + if ((XLF_FLAGV(gp) & XLF_TOFIX) != 0) { + int ud = 0; /* Update point flag */ + int th = 0; /* Try harder flag */ + sp = gp + p->oauxvv; /* Source */ + dp = gp + p->oauxv; /* Destination */ + rp = gp + p->oauxr; /* Range */ + for (f = 0; f < p->daux; f++) { + double v, ov; + v = sp[f]; + if (v < rp[2 * f]) { /* Limit new aux to valid locus range */ + if (tharder) + th = 1; + else + v = rp[2 * f]; + } + if (v > rp[2 * f + 1]) { + if (tharder) + th = 1; + else + v = rp[2 * f + 1]; + } + ov = dp[f]; /* Old aux value */ + if (fabs(ov - v) > 0.001) { + dp[f] = (float)v; + ud = 1; /* Worth updating it */ + } + } + + if (ud != 0) { + XLF_FLAGV(gp) |= XLF_UPDATE; + if (th != 0) + XLF_FLAGV(gp) |= XLF_HARDER; + } + if (XLF_FLAGV(gp) & XLF_UPDATE) + tcount++; + } + + /* Increment the counter */ + for (ce = 0; ce < p->din; ce++) { + if (++coa[ce] < p->cmax[ce]) + break; /* No carry */ + coa[ce] = p->cmin[ce]; + } + } +#ifdef DO_STATS + printf("~1 totol no. cells that will change = %d\n",tcount); +#endif + + return tcount; +} + +/* ----------------------------------------- */ +/* Update the grid values given the new auxiliary targets. */ +/* Reset the TOFIX flags. */ +static void fix_grid( +xifs *p, +double auxw /* Compromise PCS factor */ +) { + int e, f; + int coa[MAX_CHAN], ce; /* grid counter */ + + /* Init the grid counter */ + for (ce = 0; ce < p->din; ce++) + coa[ce] = p->cmin[ce]; + ce = 0; + + /* Itterate throught the PCS clut grid cells, */ + /* computing new auxiliary values */ + while (ce < p->din) { + float *gp, *ep; /* Grid pointer */ + + /* Compute base of cell pointer */ + gp = p->g; /* Grid pointer */ + for (e = 0; e < p->din; e++) + gp += coa[e] * p->fci[e]; + + if ((XLF_FLAGV(gp) & XLF_TOFIX) != 0) { + if ((XLF_FLAGV(gp) & XLF_UPDATE) != 0) { + double out[MAX_CHAN], auxv[MAX_CHAN], pcs[MAX_CHAN]; + double auxrv[MAX_CHAN]; + + ep = gp + p->opcs; + for (f = 0; f < 3; f++) /* Get PCS values */ + pcs[f] = (double)ep[f]; + + ep = gp + p->oauxv; + for (f = 0; f < p->daux; f++) /* Get auxiliary values */ + auxv[f] = (double)ep[f]; + + ep = gp + p->oout; + + if ((XLF_FLAGV(gp) & XLF_HARDER) != 0) { + /* Use "try harder" PCS->devicve lookup */ + + /* Set starting value */ + for (f = 0; f < p->dout; f++) + out[f] = (double)ep[f]; + + if (comp_pcs(p, auxw, auxrv, auxv, pcs, out) == 0) { + for (f = 0; f < p->dout; f++) /* Save the new output values */ + ep[f] = (float)out[f]; + } +#ifndef NEVER + else { + printf("~9 comp_pcs failed!\n"); + } +#endif + + /* Save actual auxiliary values */ + ep = gp + p->oauxv; + for (f = 0; f < p->daux; f++) + ep[f] = (float)auxrv[f]; + + } else { + + /* Lookup output value with different auxiliary target */ + p->clutpcsfunc(p->cbctx, out, auxv, pcs); + + for (f = 0; f < p->dout; f++) /* Save the new output values */ + ep[f] = (float)out[f]; + + /* assume auxiliary target will have been met */ + } + } + + XLF_FLAGV(gp) &= ~(XLF_TOFIX | XLF_UPDATE | XLF_HARDER); + } + + /* Increment the counter */ + for (ce = 0; ce < p->din; ce++) { + if (++coa[ce] < p->cmax[ce]) + break; /* No carry */ + coa[ce] = p->cmin[ce]; + } + } +} + +/* ----------------------------------------- */ + +/* minimizer callback function */ +static double minfunc( /* Return function value */ +void *fdata, /* Opaque data pointer */ +double *tp /* Device input value */ +) { + xifs *p = (xifs *)fdata; + double pcs[3]; + double auxv[MAX_CHAN]; + double olimit; + double fval; + double tval; + int e, f; + +#define GWHT 1000.0 + + /* Check if the device input values are outside (assumed) device range 0.0 - 1.0 */ + for (tval = 0.0, f = 0; f < p->dout; f++) { + if (tp[f] < 0.0) { + if (tp[0] > tval) + tval = tp[0]; + tp[f] = 0.0; + } else if (tp[f] > 1.0) { + if ((tp[0] - 1.0) > tval) + tval = (tp[0] - 1.0); + tp[f] = 1.0; + } + } + + /* Figure PCS, auxiliary error, and amount over ink limit */ + p->clutipcsfunc(p->cbctx, pcs, &olimit, auxv, tp); + + if (olimit > tval) + tval = olimit; + + fval = GWHT * tval; /* Largest value that exceeds device/ink range */ + + /* Figure auxiliary error */ + for (tval = 0.0, e = 0; e < p->daux; e++) { + double tt; + tt = auxv[e] - p->m_auxv[e]; + tval += tt * tt; + } + + fval += p->m_auxw * sqrt(tval); + + /* Figure PCS error */ + for (tval = 0.0, e = 0; e < 3; e++) { + double tt; + tt = pcs[e] - p->m_pcs[e]; + tval += tt * tt; + } + + fval += sqrt(tval); + +//printf("~9 minfunc returning error %f on %f %f %f %f\n", fval, tp[0], tp[1], tp[2], tp[3]); + return fval; +} + +/* Given a PCS target, and auxiliary target, and a current */ +/* device value, return a device value that is a better */ +/* tradeoff between the PCS target and the auxiliary target. */ +/* return non-zero on error */ +static int comp_pcs( +xifs *p, /* Aux fix structure */ +double auxw, /* Auxiliary error weighting factor (ie. 5 - 100) */ +double *auxrv, /* If not NULL, return actual auxiliary acheived */ +double *auxv, /* Auxiliary target value */ +double *pcs, /* PCS target value */ +double *dev /* Device starting value, return value */ +) { + int i; + double rv; + double ss[MAX_CHAN]; + double check[3]; /* Check PCS */ +#ifdef NEVER /* Diagnostic info */ +double start[3]; /* Starting PCS */ +double auxst[1]; /* Starting aux */ +double auxch[1]; /* Auxiliary check value */ +p->clutipcsfunc(p->cbctx, start, NULL, auxst, dev); +#endif + + p->m_auxw = auxw; + + for (i = 0; i < p->daux; i++) + p->m_auxv[i] = auxv[i]; + + for (i = 0; i < 3; i++) + p->m_pcs[i] = pcs[i]; + + /* Set initial search size */ + for (i = 0; i < p->dout; i++) + ss[i] = 0.3; + + if (powell(&rv, p->dout, dev, ss, 0.001, 1000, minfunc, p, NULL, NULL) != 0) { + return 1; + } + + /* Sanitise device values */ + for (i = 0; i < p->dout; i++) { + double v = dev[i]; + if (v < 0.0) + v = 0.0; + else if (v > 1.0) + v = 1.0; + dev[i] = v; + } + + if (auxrv != NULL) { /* Return actual auxiliary */ + p->clutipcsfunc(p->cbctx, check, NULL, auxrv, dev); + } + +#ifdef NEVER /* Diagnostic info */ +p->clutipcsfunc(p->cbctx, check, NULL, auxch, dev); +printf("~9 comp_pcs target aux %f PCS %f %f %f\n", auxv[0], pcs[0], pcs[1], pcs[2]); +printf("~9 returning error %f on %f %f %f %f\n", rv, dev[0], dev[1], dev[2], dev[3]); +printf("~9 PCS start = %f %f %f (%f)\n",start[0], start[1], start[2], + absdiffn(3, start, pcs)); +printf("~9 PCS finish = %f %f %f (%f)\n",check[0], check[1], check[2], + absdiffn(3, check, pcs)); +printf("~9 PCS delta = %f, aux delta = %f\n", absdiffn(3, start, check), + fabs(auxst[0] - auxch[0])); +#endif /* NEVER */ + + return 0; +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + |