From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- spectro/dispcal.c | 5356 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 5356 insertions(+) create mode 100644 spectro/dispcal.c (limited to 'spectro/dispcal.c') diff --git a/spectro/dispcal.c b/spectro/dispcal.c new file mode 100644 index 0000000..accb79a --- /dev/null +++ b/spectro/dispcal.c @@ -0,0 +1,5356 @@ + +/* + * Argyll Color Correction System + * Display callibrator. + * + * Author: Graeme W. Gill + * Date: 14/10/2005 + * + * Copyright 1996 - 2013 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 program displays test patches, and takes readings from a display device */ +/* in order to create a RAMDAC calibration curve (usually stored in the ICC vcgt tag) */ + +/* This is the third version of the program. */ + +/* TTBD + + Try to improve calibration speed by using adaptive + measurement set, rather than fixed resolution doubling ? + (ie. just measure at troublesome points using a "divide in half" + strategy ?. Estimate error between measurement points and + pick the next largest error.) + + Add option to use L*u*v* DE's, as this is used in + some video standards. They sometime use u*v* as + a color tollerance too (see EBU TECH 3320). + + Add a white point option that makes the target the + closest temperature to the native one of the display :- + ie. it moves the display to the closest point on the + chosen locus to RGB 1,1,1. + ie. should it do this if "-t" or "-T" + with no specific temperature is chosen ? + + Change white point gamut clipping to be a measurement + search rather than computing from primary XYZ ? + + Add bell at end of calibration ? + + + Add option to plot graph of native and calibrated RGB ? + + Add a "delta E" number to the interactive adjustments, + so the significance of the error can be judged ? + + Need to add flare measure/subtract, to improve + projector calibration ? - need to add to dispread too. + + Instead of measuring/matching output at 50% device input as + measure of gamma, what about inverting it - measure/match device + values at 50% perceptual (18%) output value ? + [ Hmm. Current method is OK because a good perceptual + display gives about 18% output at 50% device input.] + + + The verify (-E) may not be being done correctly. + Like update, shouldn't it read the .cal file to set what's + being calibrated aganist ? (This would fix missing ambient value too!) + + What about the "Read the base test set" - aren't + there numbers then used to tweak the black aim point + in "Figure out the black point target" - Yes they are !! + Verify probably shouldn't work this way. + + */ + +#ifdef __MINGW32__ +# define WINVER 0x0500 +#endif + +#include +#include +#include +#include +#include +#include +#include +#if defined (NT) +#include +#endif +#include "copyright.h" +#include "aconfig.h" +#include "numlib.h" +#include "xicc.h" +#include "xspect.h" +#include "xcolorants.h" +#include "ccmx.h" +#include "ccss.h" +#include "cgats.h" +#include "insttypes.h" +#include "conv.h" +#include "icoms.h" +#include "inst.h" +#include "dispwin.h" +#include "dispsup.h" +#include "rspl.h" +#include "moncurve.h" +#include "targen.h" +#include "ofps.h" +#include "icc.h" +#include "sort.h" +#include "instappsup.h" +#include "spyd2setup.h" /* Enable Spyder 2 access */ + +#undef DEBUG +#undef DEBUG_OFFSET /* Keep test window out of the way */ +#undef DEBUG_PLOT /* Plot curve each time around */ +#undef CHECK_MODEL /* Do readings to check the accuracy of our model */ +#undef SHOW_WINDOW_ONFAKE /* Display a test window up for a fake device */ + +/* Invoke with -dfake for testing with a fake device. */ +/* Will use a fake.icm/.icc profile if present, or a built in fake */ +/* device behaviour if not. */ + +#define COMPORT 1 /* Default com port 1..4 */ +#define OPTIMIZE_MODEL /* Adjust model for best fit */ +#define REFINE_GAIN 0.80 /* Refinement correction damping/gain */ +#define MAX_RPTS 12 /* Maximum tries at making a sample meet the current threshold */ +#define VER_RES 100 /* Verification resolution */ +#define NEUTRAL_BLEND_RATE 4.0 /* Default rate of transition for -k factor < 1.0 (power) */ +#define ADJ_JACOBIAN /* Adjust the Jacobian predictor matrix each time */ +#define JAC_COR_FACT 0.4 /* Amount to correct Jacobian by (to filter noise) */ +#define REMEAS_JACOBIAN /* Re-measure Jacobian */ +#define MOD_DIST_POW 1.6 /* Power used to distribute test samples for model building */ +#define REFN_DIST_POW 1.6 /* Power used to distribute test samples for grey axis refinement */ +#define CHECK_DIST_POW 1.6 /* Power used to distribute test samples for grey axis checking */ +#define THRESH_SCALE_POW 0.5 /* Amount to loosen threshold for first itterations */ +#define CAL_RES 256 /* Resolution of calibration table to produce. */ +#define CLIP /* Clip RGB during refinement */ +#define RDAC_SMOOTH 0.3 /* RAMDAC curve fitting smoothness */ +#define MEAS_RES /* Measure the RAMNDAC entry size */ + +#ifdef DEBUG_PLOT +#include "plot.h" +#endif + +#if defined(DEBUG) + +#define DBG(xxx) fprintf xxx ; +#define dbgo stderr +#else +#define DBG(xxx) +#endif /* DEBUG */ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Sample points used in initial device model optimisation */ + +typedef struct { + double dev[3]; /* Device values */ + double lab[3]; /* Read value */ + double w; /* Weighting */ +} optref; + +/* - - - - - - - - - - - - - - - - - - - */ +/* device RGB inverse solution code */ + +/* Selected transfer curve */ +typedef enum { + gt_power = 0, /* A simple power */ + gt_Lab = 1, /* The L* curve */ + gt_sRGB = 2, /* The sRGB curve */ + gt_Rec709 = 3, /* REC 709 video standard */ + gt_SMPTE240M = 4 /* SMTPE 240M video standard */ +} gammatype; + +/* Context for calibration solution */ +typedef struct { + double wh[3]; /* White absolute XYZ value */ + double bk[3]; /* Black absolute XYZ value */ + + /* Target model */ + gammatype gammat; /* Transfer curve type */ + double egamma; /* Effective Gamma target */ + double oofff; /* proportion of output offset vs input offset (default 1.0) */ + double gioff; /* Gamma curve input zero offset */ + double gooff; /* Target output offset (normalised to Y max of 1.0) */ + int nat; /* Flag - nz if native white target */ + double nbrate; /* Neutral blend weight (power) */ + + /* Viewing conditions adjustment */ + int vc; /* Flag, nz to enable viewing conditions adjustment */ + icxcam *svc; /* Source viewing conditions */ + icxcam *dvc; /* Destination viewing conditions */ + double vn0, vn1; /* Normalisation values */ + + double nwh[3]; /* Target white normalised XYZ value (Y = 1.0) */ + double twh[3]; /* Target white absolute XYZ value */ + icmXYZNumber twN; /* Same as above as XYZNumber */ + + double tbk[3]; /* Target black point color */ + icmXYZNumber tbN; /* Same as above as XYZNumber */ + + /* Device model */ + double fm[3][3]; /* Forward, aprox. linear RGB -> XYZ */ + double bm[3][3]; /* Backwards, aprox. XYZ -> linear RGB */ + mcv *dcvs[3]; /* Device RGB channel to linearised RGB curves */ + /* These are always normalized to map 1.0 to 1.0 */ + + /* Current state */ + mcv *rdac[3]; /* Current RGB to RGB ramdac curves */ + + double xyz[3]; /* Target xyz value */ + + /* optimisation information */ + int np; /* Total number of optimisation parameters */ + int co[3]; /* Offset in the parameters to each curve offset */ + int nc[3]; /* Number of for each curve */ + int nrp; /* Total number of reference points */ + optref *rp; /* reference points */ + double *dtin_iv; /* Temporary array :- dp for input curves */ +} calx; + +/* - - - - - - - - - - - - - - - - - - - */ +/* Ideal target curve definitions */ + +/* Convert ideal device (0..1) to target Y value (0..1) */ +static double dev2Y(calx *x, double egamma, double vv) { + + switch(x->gammat) { + case gt_power: { + vv = pow(vv, egamma); + break; + } + case gt_Lab: { + vv = icmL2Y(vv * 100.0); + break; + } + case gt_sRGB: { + if (vv <= 0.03928) + vv = vv/12.92; + else + vv = pow((0.055 + vv)/1.055, 2.4); + break; + } + case gt_Rec709: { + if (vv <= 0.081) + vv = vv/4.5; + else + vv = pow((0.099 + vv)/1.099, 1.0/0.45); + break; + } + case gt_SMPTE240M: { + if (vv <= 0.0913) + vv = vv/4.0; + else + vv = pow((0.1115 + vv)/1.1115, 1.0/0.45); + break; + } + default: + error("Unknown gamma type"); + } + return vv; +} + +/* Convert target Y value (0..1) to ideal device (0..1) */ +static double Y2dev(calx *x, double egamma, double vv) { + + switch(x->gammat) { + case gt_power: { + vv = pow(vv, 1.0/egamma); + break; + } + case gt_Lab: { + vv = icmY2L(vv) * 0.01; + break; + } + case gt_sRGB: { + if (vv <= 0.00304) + vv = vv * 12.92; + else + vv = pow(vv, 1.0/2.4) * 1.055 - 0.055; + break; + } + case gt_Rec709: { + if (vv <= 0.018) + vv = vv * 4.5; + else + vv = pow(vv, 0.45) * 1.099 - 0.099; + break; + } + case gt_SMPTE240M: { + if (vv <= 0.0228) + vv = vv * 4.0; + else + vv = pow(vv, 0.45) * 1.1115 - 0.1115; + break; + } + default: + error("Unknown gamma type"); + } + return vv; +} + +/* - - - - - - - - - - - - - - - - - - - */ +/* Compute a viewing environment Y transform */ + +static double view_xform(calx *x, double in) { + double out = in; + + if (x->vc != 0) { + double xyz[3], Jab[3]; + + xyz[0] = in * x->nwh[0]; /* Compute value on neutral axis */ + xyz[1] = in * x->nwh[1]; + xyz[2] = in * x->nwh[2]; + x->svc->XYZ_to_cam(x->svc, Jab, xyz); + x->dvc->cam_to_XYZ(x->dvc, xyz, Jab); + + out = xyz[1] * x->vn1 + x->vn0; /* Apply scaling factors */ + } + return out; +} + +/* - - - - - - - - - - - - - - - - - - - */ + +/* Info for optimization */ +typedef struct { + double thyr; /* 50% input target */ + double roo; /* 0% input target */ +} gam_fits; + +/* gamma + input offset function handed to powell() */ +static double gam_fit(void *dd, double *v) { + gam_fits *gf = (gam_fits *)dd; + double gamma = v[0]; + double ioff = v[1]; + double rv = 0.0; + double tt; + + if (gamma < 0.0) { + rv += 100.0 * -gamma; + gamma = 0.0; + } + if (ioff < 0.0) { + rv += 100.0 * -ioff; + ioff = 0.0; + } else if (ioff > 0.999) { + rv += 100.0 * (ioff - 0.999); + ioff = 0.999; + } + tt = gf->roo - pow(ioff, gamma); + rv += tt * tt; + tt = gf->thyr - pow(0.5 + (1.0 - 0.5) * ioff, gamma); + rv += tt * tt; + +//printf("~1 gam_fit %f %f returning %f\n",ioff,gamma,rv); + return rv; +} + + +/* Given the advertised gamma and the output offset, compute the */ +/* effective gamma and input offset needed. */ +/* Return the expected output value for 50% input. */ +/* (It's assumed that gooff is normalised the target brightness) */ +static double tech_gamma( + calx *x, + double *pegamma, /* return effective gamma needed */ + double *pooff, /* return output offset needed */ + double *pioff, /* return input offset needed */ + double egamma, /* effective gamma needed (> 0.0 if valid, overrides gamma) */ + double gamma, /* advertised gamma needed */ + double tooff /* Total ouput offset needed */ +) { + int i; + double rv; + double gooff = 0.0; /* The output offset applied */ + double gioff = 0.0; /* The input offset applied */ + double roo; /* Remaining output offset accounted for by input offset */ + + /* Compute the output offset that will be applied */ + gooff = tooff * x->oofff; + roo = (tooff - gooff)/(1.0 - gooff); + +//printf("~1 gooff = %f, roo = %f\n",gooff,roo); + + /* Now compute the input offset that will be needed */ + if (x->gammat == gt_power && egamma <= 0.0) { + gam_fits gf; + double op[2], sa[2], rv; + + gf.thyr = pow(0.5, gamma); /* Advetised 50% target */ + gf.thyr = (gf.thyr - gooff)/(1.0 - gooff); /* Target before gooff is added */ + gf.roo = roo; + + op[0] = gamma; + op[1] = pow(roo, 1.0/gamma); + sa[0] = 0.1; + sa[1] = 0.01; + + if (powell(&rv, 2, op, sa, 1e-6, 500, gam_fit, (void *)&gf, NULL, NULL) != 0) + warning("Computing effective gamma and input offset is inaccurate"); + + if (rv > 1e-5) { + warning("Computing effective gamma and input offset is inaccurate (%f)",rv); + } + egamma = op[0]; + gioff = op[1]; + +//printf("~1 Result gioff %f, gooff %f, egamma %f\n",gioff, gooff, egamma); +//printf("~1 Verify 0.0 in -> out = %f, tooff = %f\n",gooff + dev2Y(x, egamma, gioff) * (1.0 - gooff),tooff); +//printf("~1 Verify 0.5 out = %f, target %f\n",gooff + dev2Y(x, egamma, gioff + 0.5 * (1.0 - gioff)) * (1.0 - gooff), pow(0.5, gamma)); + + } else { + gioff = Y2dev(x, egamma, roo); +//printf("~1 Result gioff %f, gooff %f\n",gioff, gooff); +//printf("~1 Verify 0.0 in -> out = %f, tooff = %f\n",gooff + dev2Y(x, egamma, gioff) * (1.0 - gooff),tooff); + } + + /* Compute the 50% output value */ + rv = gooff + dev2Y(x, egamma, gioff + 0.5 * (1.0 - gioff)) * (1.0 - gooff); + + if (pegamma != NULL) + *pegamma = egamma; + if (pooff != NULL) + *pooff = gooff; + if (pioff != NULL) + *pioff = gioff; + return rv; +} + +/* Compute approximate advertised gamma from black/50% grey/white readings, */ +/* (assumes a zero based gamma curve shape) */ +static double pop_gamma(double bY, double gY, double wY) { + int i; + double grat, brat, gioff, gvv, gamma; + + grat = gY/wY; + brat = bY/wY; + + gamma = log(grat) / log(0.5); + return gamma; +} + +/* - - - - - - - - - - - - - - - - - - - */ + +/* Return the xyz that is predicted by our aproximate device model */ +/* by the given device RGB. */ +static void fwddev(calx *x, double xyz[3], double rgb[3]) { + double lrgb[3]; + int j; + +//printf("~1 fwddev called with rgb %f %f %f\n",rgb[0],rgb[1],rgb[2]); + + /* Convert device RGB into linear light RGB via curves */ + for (j = 0; j < 3; j++) + lrgb[j] = x->dcvs[j]->interp(x->dcvs[j], rgb[j]); + +//printf("~1 fwddev got linear RGB %f %f %f\n",lrgb[0],lrgb[1],lrgb[2]); + + /* Convert linear light RGB into XYZ via the matrix */ + icmMulBy3x3(xyz, x->fm, lrgb); + +//printf("~1 fwddev got final xyz %f %f %f\n",xyz[0],xyz[1],xyz[2]); +} + +/* Return the closest device RGB predicted by our aprox. device model */ +/* to generate the given xyz. */ +static void invdev(calx *x, double rgb[3], double xyz[3]) { + double lrgb[3]; + int j; + +//printf("~1 invdev called with xyz %f %f %f\n",xyz[0],xyz[1],xyz[2]); + + /* Convert XYZ to linear light RGB via the inverse matrix */ + icmMulBy3x3(lrgb, x->bm, xyz); +//printf("~1 invdev; lin light rgb = %f %f %f\n",lrgb[0],lrgb[1],lrgb[2]); + + /* Convert linear light RGB to device RGB via inverse curves */ + for (j = 0; j < 3; j++) { + lrgb[j] = x->dcvs[j]->inv_interp(x->dcvs[j], lrgb[j]); + if (lrgb[j] < 0.0) { +#ifdef CLIP + lrgb[j] = 0.0; +#endif + } else if (lrgb[j] > 1.0) { +#ifdef CLIP + lrgb[j] = 1.0; +#endif + } + } +//printf("~1 invdev; inverse curves rgb = %f %f %f\n",lrgb[0],lrgb[1],lrgb[2]); + if (rgb != NULL) { + rgb[0] = lrgb[0]; + rgb[1] = lrgb[1]; + rgb[2] = lrgb[2]; + } +} + +/* Return the closest linear device RGB predicted by our aprox. device matrix */ +/* to generate the given xyz. */ +/* Return > 0 if clipped */ +static double invlindev(calx *x, double rgb[3], double xyz[3]) { + double lrgb[3]; + double clip = 0.0; + int j; + +//printf("~1 invlindev called with xyz %f %f %f\n",xyz[0],xyz[1],xyz[2]); + + /* Convert XYZ to linear light RGB via the inverse matrix */ + icmMulBy3x3(lrgb, x->bm, xyz); +//printf("~1 invlindev; lin light rgb = %f %f %f\n",lrgb[0],lrgb[1],lrgb[2]); + + /* Check for out of gamut */ + for (j = 0; j < 3; j++) { + if (lrgb[j] < 0.0) { + if (-lrgb[j] > clip) + clip = -lrgb[j]; + lrgb[j] = 0.0; + } else if (lrgb[j] > 1.0) { + if ((lrgb[j]-1.0) > clip) + clip = (lrgb[j]-1.0); + lrgb[j] = 1.0; + } + } +//printf("~1 invlindev; clipped rgb = %f %f %f, clip = %f \n",lrgb[0],lrgb[1],lrgb[2],clip); + if (rgb != NULL) { + rgb[0] = lrgb[0]; + rgb[1] = lrgb[1]; + rgb[2] = lrgb[2]; + } + return clip; +} + +/* Overall optimisation support */ + +/* Set the optimsation parameter number and offset values in calx, */ +/* and return an array filled in with the current parameters. */ +/* Allocate temporary arrays */ +static double *dev_get_params(calx *x) { + double *p, *tp; + int i, j; + + x->np = 9; + for (i = 0; i < 3; i++) + x->np += x->dcvs[i]->luord; + + if ((p = (double *)malloc(x->np * sizeof(double))) == NULL) + error("dev_params malloc failed"); + + tp = p; + + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + *tp++ = x->fm[i][j]; + + for (i = 0; i < 3; i++) { + x->co[i] = tp - p; /* Offset to start */ + for (j = 0; j < x->dcvs[i]->luord; j++) + *tp++ = x->dcvs[i]->pms[j]; + x->nc[i] = (tp - p) - x->co[i]; /* Number */ + } + + if ((x->dtin_iv = (double *)malloc(x->np * sizeof(double))) == NULL) + error("dev_params malloc failed"); + + return p; +} + +/* Given a set of parameters, put them back into the model */ +/* Normalize them so that the curve maximum is 1.0 too. */ +static void dev_put_params(calx *x, double *p) { + int i, j; + double scale[3]; + + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + x->fm[i][j] = *p++; + for (i = 0; i < 3; i++) + for (j = 0; j < x->dcvs[i]->luord; j++) + x->dcvs[i]->pms[j] = *p++; + + /* Figure out how we have to scale the curves */ + for (j = 0; j < 3; j++) { + scale[j] = x->dcvs[j]->interp(x->dcvs[j], 1.0); + x->dcvs[j]->force_scale(x->dcvs[j], 1.0); + } + + /* Scale the matrix to compensate */ + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + x->fm[i][j] *= scale[j]; +} + +/* Device model optimisation function handed to powell() */ +static double dev_opt_func(void *edata, double *v) { + calx *x = (calx *)edata; + int i, j; + double tw = 0.0; + double rv, smv; + +#ifdef NEVER + printf("params ="); + for (i = 0; i < x->np; i++) + printf(" %f",v[i]); + printf("\n"); +#endif + + /* For all our data points */ + rv = 0.0; + for (i = 0; i < x->nrp; i++) { + double lrgb[3]; /* Linear light RGB */ + double xyz[3], lab[3]; + double de; + + /* Convert through device curves */ + for (j = 0; j < 3; j++) + lrgb[j] = x->dcvs[j]->interp_p(x->dcvs[j], v + x->co[j], x->rp[i].dev[j]); + + /* Convert linear light RGB into XYZ via the matrix */ + icxMulBy3x3Parm(xyz, v, lrgb); + + /* Convert to Lab */ + icmXYZ2Lab(&x->twN, lab, xyz); + + /* Compute delta E squared */ + de = icmCIE94sq(lab, x->rp[i].lab) * x->rp[i].w; +#ifdef NEVER + printf("point %d DE %f, Lab is %f %f %f, should be %f %f %f\n", +i, sqrt(de), lab[0], lab[1], lab[2], x->rp[i].lab[0], x->rp[i].lab[1], x->rp[i].lab[2]); +#endif + rv += de; + tw += x->rp[i].w; + } + + /* Normalise error to be a weighted average delta E squared and scale smoothing */ + rv /= (tw * 5.0); + + /* Sum with shaper parameters squared, to */ + /* minimise unsconstrained "wiggles" */ + smv = 0.0; + for (j = 0; j < 3; j++) + smv += x->dcvs[j]->shweight_p(x->dcvs[j], v + x->co[j], 1.0); + rv += smv; + +#ifdef NEVER + printf("rv = %f (%f)\n",rv, smv); +#endif + return rv; +} + + +/* Device model optimisation function handed to conjgrad() */ +static double dev_dopt_func(void *edata, double *dv, double *v) { + calx *x = (calx *)edata; + int i, j, k; + int f, ee, ff, jj; + double tw = 0.0; + double rv, smv; + + double dmato_mv[3][9]; /* Del in mat out due to del in matrix param vals */ + double dmato_tin[3][3]; /* Del in mat out due to del in matrix input values */ + double dout_lab[3][3]; /* Del in out due to XYZ to Lab conversion */ + double de_dout[2][3]; /* Del in delta E due to input Lab values */ + +#ifdef NEVER + printf("params ="); + for (i = 0; i < x->np; i++) + printf(" %f",v[i]); + printf("\n"); +#endif + + /* Zero the accumulated partial derivatives */ + for (i = 0; i < x->np; i++) + dv[i] = 0.0; + + /* For all our data points */ + rv = 0.0; + for (i = 0; i < x->nrp; i++) { + double lrgb[3]; /* Linear light RGB */ + double xyz[3], lab[3]; + + /* Apply the input channel curves */ + for (j = 0; j < 3; j++) + lrgb[j] = x->dcvs[j]->dinterp_p(x->dcvs[j], v + x->co[j], + x->dtin_iv + x->co[j], x->rp[i].dev[j]); + + /* Convert linear light RGB into XYZ via the matrix */ + icxdpdiMulBy3x3Parm(xyz, dmato_mv, dmato_tin, v, lrgb); + + /* Convert to Lab */ + icxdXYZ2Lab(&x->twN, lab, dout_lab, xyz); + + /* Compute delta E squared */ +//printf("~1 point %d: Lab is %f %f %f, should be %f %f %f\n", +//i, lab[0], lab[1], lab[2], x->rp[i].lab[0], x->rp[i].lab[1], x->rp[i].lab[2]); + rv += icxdCIE94sq(de_dout, lab, x->rp[i].lab) * x->rp[i].w; + de_dout[0][0] *= x->rp[i].w; + de_dout[0][1] *= x->rp[i].w; + de_dout[0][2] *= x->rp[i].w; + tw += x->rp[i].w; + + /* Compute and accumulate partial difference values for each parameter value */ + + /* Input channel curves */ + for (ee = 0; ee < 3; ee++) { /* Parameter input chanel */ + for (k = 0; k < x->nc[ee]; k++) { /* Param within channel */ + double vv = 0.0; + jj = x->co[ee] + k; /* Overall input curve param */ + + for (ff = 0; ff < 3; ff++) { /* Lab channels */ + for (f = 0; f < 3; f++) { /* XYZ channels */ + vv += de_dout[0][ff] * dout_lab[ff][f] + * dmato_tin[f][ee] * x->dtin_iv[jj]; + } + } + dv[jj] += vv; + } + } + + /* Matrix parameters */ + for (k = 0; k < 9; k++) { /* Matrix parameter */ + double vv = 0.0; + + for (ff = 0; ff < 3; ff++) { /* Lab channels */ + for (f = 0; f < 3; f++) { /* XYZ channels */ + vv += de_dout[0][ff] * dout_lab[ff][f] + * dmato_mv[f][k]; + } + } + dv[k] += vv; + } + } + + /* Normalise error to be a weighted average delta E squared and scale smoothing */ + rv /= (tw * 1200.0); + for (i = 0; i < x->np; i++) + dv[i] /= (tw * 900.0); + + /* Sum with shaper parameters squared, to */ + /* minimise unsconstrained "wiggles" */ + smv = 0.0; + for (j = 0; j < 3; j++) + smv += x->dcvs[j]->dshweight_p(x->dcvs[j], v + x->co[j], x->dtin_iv + x->co[j], 1.0); + rv += smv; + +#ifdef NEVER + printf("drv = %f (%f)\n",rv, smv); +#endif + return rv; +} + +#ifdef NEVER +/* Check partial derivative function within dev_opt_func() using powell() */ + +static double dev_opt_func(void *edata, double *v) { + calx *x = (calx *)edata; + int i; + double dv[2000]; + double rv, drv; + double trv; + + rv = dev_opt_func_(edata, v); + drv = dev_dopt_func(edata, dv, v); + + if (fabs(rv - drv) > 1e-6) { + printf("######## RV MISMATCH is %f should be %f ########\n",drv, rv); + exit(0); + } + + /* Check each parameter delta */ + for (i = 0; i < x->np; i++) { + double del; + + v[i] += 1e-7; + trv = dev_opt_func_(edata, v); + v[i] -= 1e-7; + + /* Check that del is correct */ + del = (trv - rv)/1e-7; + if (fabs(dv[i] - del) > 1.0) { +//printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv); + printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del); + exit(0); + } + } + return rv; +} +#endif + +/* =================================================================== */ + +/* White point brightness optimization function handed to powell. */ +/* Maximize brigtness while staying within gamut */ +static double wp_opt_func(void *edata, double *v) { + calx *x = (calx *)edata; + double wxyz[3], rgb[3]; + int j; + double rv = 0.0; + + wxyz[0] = v[0] * x->twh[0]; + wxyz[1] = v[0] * x->twh[1]; + wxyz[2] = v[0] * x->twh[2]; + +//printf("~1 wp_opt_func got scale %f, xyz = %f %f %f\n", +//v[0],wxyz[0],wxyz[1],wxyz[2]); + + if ((rv = invlindev(x, rgb, wxyz)) > 0.0) { /* Out of gamut */ + rv *= 1e5; +//printf("~1 out of gamut %f %f %f returning %f\n", rgb[0], rgb[1], rgb[2], rv); + return rv; + } + /* Maximize scale factor */ + if (v[0] < 0.00001) + rv = 1.0/0.00001; + else + rv = 1.0/v[0]; + +//printf("~1 %f %f %f returning %f\n", rgb[0], rgb[1], rgb[2], rv); + return rv; +} + +/* =================================================================== */ +/* Structure to save aproximate model readings in */ +typedef struct { + double v; /* Input value */ + double xyz[3]; /* Reading */ +} sxyz; + +/* ------------------------------------------------------------------- */ +#if defined(__APPLE__) && defined(__POWERPC__) + +/* Workaround for a ppc gcc 3.3 optimiser bug... */ +/* It seems to cause a segmentation fault instead of */ +/* converting an integer loop index into a float, */ +/* when there are sufficient variables in play. */ +static int gcc_bug_fix(int i) { + static int nn; + nn += i; + return nn; +} +#endif /* APPLE */ + + +/* =================================================================== */ +/* Calibration sample point support. This allows the successive */ +/* refinement of our neutral sample points */ + +/* A sample point */ +typedef struct { + double v; /* Desired input value */ + double rgb[3]; /* Input value through calibration curves */ + double tXYZ[3]; /* Target XYZ */ + double XYZ[3]; /* Read XYZ */ + double deXYZ[3]; /* Delta XYZ wanted to target */ + double de; /* Delta Lab to neutral target */ + double dc; /* Delta XYZ to neutral target */ + double peqde; /* Delta Lab to previous point value (last pass) */ + double hde; /* Hybrid de composed of de and peqde */ + + double pXYZ[3]; /* Previous measured XYZ */ + double pdXYZ[3]; /* Delta XYZ intended from previous measure */ + double pdrgb[3]; /* Delta rgb made to previous */ + + double dXYZ[3]; /* Actual delta XYZ resulting from previous delta rgb */ + + double j[3][3]; /* Aproximate Jacobian (del RGB -> XYZ) */ + double ij[3][3]; /* Aproximate inverse Jacobian (del XYZ-> del RGB) */ + double fb_ij[3][3]; /* Copy of initial inverse Jacobian, used as a fallback */ +} csp; + +/* All the sample points */ +typedef struct { + int no; /* Number of samples */ + int _no; /* Allocation */ + csp *s; /* List of samples */ +} csamp; + +static void free_alloc_csamp(csamp *p) { + if (p->s != NULL) + free(p->s); + p->s = NULL; +} + +/* Initialise v values */ +static void init_csamp_v(csamp *p, calx *x, int psrand) { + int i, j; + sobol *so = NULL; + + if (psrand != 0) { /* Use pseudo random distribution for verification */ + if ((so = new_sobol(1)) == NULL) + error("New sobol failed"); + } + + /* Generate the sample points */ + for (i = 0; i < p->no; i++) { + double vv; + +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(i); +#endif + if (so != NULL) { + if (i == 0) + vv = 1.0; + else if (i == 1) + vv = 0.0; + else + so->next(so, &vv); + } else + vv = i/(p->no - 1.0); + vv = pow(vv, REFN_DIST_POW); /* Skew sample points to be slightly perceptual */ + p->s[i].v = vv; + } + + if (so != NULL) { + /* Sort it so white is last */ +#define HEAP_COMPARE(A,B) (A.v < B.v) + HEAPSORT(csp,p->s,p->no) +#undef HEAP_COMPARE + so->del(so); + } +} + +/* Initialise txyz values from v values */ +static void init_csamp_txyz(csamp *p, calx *x, int fixdev) { + int i, j; + double tbL[3]; /* tbk as Lab */ + + /* Convert target black from XYZ to Lab here, */ + /* in case twN has changed at some point. */ + icmXYZ2Lab(&x->twN, tbL, x->tbk); + + /* Set the sample points targets */ + for (i = 0; i < p->no; i++) { + double y, vv; + double XYZ[3]; /* Existing XYZ value */ + double Lab[3]; + double bl; + + vv = p->s[i].v; + + /* Compute target relative Y value for this device input. */ + /* We allow for any input and/or output offset */ + y = x->gooff + dev2Y(x, x->egamma, x->gioff + vv * (1.0 - x->gioff)) * (1.0 - x->gooff); + + /* Add viewing environment transform */ + y = view_xform(x, y); + + /* Convert Y to L* */ + Lab[0] = icmY2L(y); + Lab[1] = Lab[2] = 0.0; /* Target is neutral */ + + /* Compute blended neutral target a* b* */ + bl = pow((1.0 - vv), x->nbrate); /* Crossover near the black */ + Lab[1] = bl * tbL[1]; + Lab[2] = bl * tbL[2]; + + icmAry2Ary(XYZ, p->s[i].tXYZ); /* Save the existing values */ + icmLab2XYZ(&x->twN, p->s[i].tXYZ, Lab); /* New XYZ Value to aim for */ + +#ifdef DEBUG + printf("%d: target XYZ %.2f %.2f %.2f, Lab %.2f %.2f %.2f\n",i, p->s[i].tXYZ[0],p->s[i].tXYZ[1],p->s[i].tXYZ[2], Lab[0],Lab[1],Lab[2]); +#endif + } +} + + +/* Allocate the sample points and initialise them with the */ +/* target device and XYZ values, and first cut device values. */ +static void init_csamp(csamp *p, calx *x, int doupdate, int verify, int psrand, int no) { + int i, j; + + p->_no = p->no = no; + + if ((p->s = (csp *)malloc(p->_no * sizeof(csp))) == NULL) + error("csamp malloc failed"); + + /* Compute v and txyz */ + init_csamp_v(p, x, psrand); + init_csamp_txyz(p, x, 0); + + /* Generate the sample points */ + for (i = 0; i < no; i++) { + double dd, vv; + +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(i); +#endif + vv = p->s[i].v; + + if (verify == 2) { /* Verifying through installed curve */ + /* Make RGB values the input value */ + p->s[i].rgb[0] = p->s[i].rgb[1] = p->s[i].rgb[2] = vv; + + } else if (doupdate) { /* Start or verify through current cal curves */ + for (j = 0; j < 3; j++) { + p->s[i].rgb[j] = x->rdac[j]->interp(x->rdac[j], vv); +#ifdef CLIP + if (p->s[i].rgb[j] < 0.0) + p->s[i].rgb[j] = 0.0; + else if (p->s[i].rgb[j] > 1.0) + p->s[i].rgb[j] = 1.0; +#endif + } + } else { /* we have model */ + /* Lookup an initial device RGB for that target by inverting */ + /* the approximate forward device model */ + p->s[i].rgb[0] = p->s[i].rgb[1] = p->s[i].rgb[2] = vv; + invdev(x, p->s[i].rgb, p->s[i].tXYZ); + } + /* Force white to be native if native flag set */ + if (x->nat && i == (no-1)) { +//printf("~1 Forcing white rgb to be 1,1,1\n"); + p->s[i].rgb[0] = p->s[i].rgb[1] = p->s[i].rgb[2] = 1.0; + } + +//printf("~1 Inital point %d rgb %f %f %f\n",i,p->s[i].rgb[0],p->s[i].rgb[1],p->s[i].rgb[2]); + + /* Compute the approximate inverse Jacobian at this point */ + /* by taking the partial derivatives wrt to each device */ + /* channel of our aproximate forward model */ + if (verify != 2) { + double refXYZ[3], delXYZ[3]; + fwddev(x, refXYZ, p->s[i].rgb); + if (vv < 0.5) + dd = 0.02; + else + dd = -0.02; + /* Matrix organization is J[XYZ][RGB] for del RGB->del XYZ*/ + for (j = 0; j < 3; j++) { + p->s[i].rgb[j] += dd; + fwddev(x, delXYZ, p->s[i].rgb); + p->s[i].j[0][j] = (delXYZ[0] - refXYZ[0]) / dd; + p->s[i].j[1][j] = (delXYZ[1] - refXYZ[1]) / dd; + p->s[i].j[2][j] = (delXYZ[2] - refXYZ[2]) / dd; + p->s[i].rgb[j] -= dd; + } + if (icmInverse3x3(p->s[i].ij, p->s[i].j)) { + error("dispcal: inverting Jacobian failed (1)"); + } + /* Make a copy of this Jacobian in case we get an invert failure later */ + icmCpy3x3(p->s[i].fb_ij, p->s[i].ij); + } + } +} + +/* Return a linear XYZ interpolation */ +static void csamp_interp(csamp *p, double xyz[3], double v) { + int i, j; + double b; + + if (p->no < 2) + error("Calling csamp_interp with less than two existing samples"); + + /* Locate the pair surrounding our input value */ + for (i = 0; i < (p->no-1); i++) { + if (v >= p->s[i].v && v <= p->s[i+1].v) + break; + } + if (i >= (p->no-1)) + error("csamp_interp out of range"); + + b = (v - p->s[i].v)/(p->s[i+1].v - p->s[i].v); + + for (j = 0; j < 3; j++) { + xyz[j] = b * p->s[i+1].XYZ[j] + (1.0 - b) * p->s[i].XYZ[j]; + } +} + +/* Re-initialise a CSP with a new number of points. */ +/* Interpolate the device values and jacobian. */ +/* Set the current rgb from the current RAMDAC curves if not verifying */ +static void reinit_csamp(csamp *p, calx *x, int verify, int psrand, int no) { + csp *os; /* Old list of samples */ + int ono; /* Old number of samples */ + int i, j, k, m; + + if (no == p->no) + return; /* Nothing has changed */ + + os = p->s; /* Save the existing per point information */ + ono = p->no; + + init_csamp(p, x, 0, 2, psrand, no); + + p->_no = p->no = no; + + /* Interpolate the current device values */ + for (i = 0; i < no; i++) { + double vv, b; + + vv = p->s[i].v; + + /* Locate the pair surrounding our target value */ + for (j = 0; j < ono-1; j++) { + if (vv >= os[j].v && vv <= os[j+1].v) + break; + } + if (j >= (ono-1)) + error("csamp interp. out of range"); + + b = (vv - os[j].v)/(os[j+1].v - os[j].v); + + for (k = 0; k < 3; k++) { + if (verify == 2) { + + p->s[i].rgb[k] = b * os[j+1].rgb[k] + (1.0 - b) * os[j].rgb[k]; + + } else { /* Lookup rgb from current calibration curves */ + for (m = 0; m < 3; m++) { + p->s[i].rgb[m] = x->rdac[m]->interp(x->rdac[m], vv); +#ifdef CLIP + if (p->s[i].rgb[m] < 0.0) + p->s[i].rgb[m] = 0.0; + else if (p->s[i].rgb[m] > 1.0) + p->s[i].rgb[m] = 1.0; +#endif + } + } + p->s[i].XYZ[k] = b * os[j+1].XYZ[k] + (1.0 - b) * os[j].XYZ[k]; + p->s[i].deXYZ[k] = b * os[j+1].deXYZ[k] + (1.0 - b) * os[j].deXYZ[k]; + p->s[i].pXYZ[k] = b * os[j+1].pXYZ[k] + (1.0 - b) * os[j].pXYZ[k]; + p->s[i].pdrgb[k] = b * os[j+1].pdrgb[k] + (1.0 - b) * os[j].pdrgb[k]; + p->s[i].dXYZ[k] = b * os[j+1].dXYZ[k] + (1.0 - b) * os[j].dXYZ[k]; +#ifdef INTERP_JAC + for (m = 0; m < 3; m++) + p->s[i].j[k][m] = b * os[j+1].j[k][m] + (1.0 - b) * os[j].j[k][m]; +#endif + + } +#ifndef INTERP_JAC + /* Create a Jacobian at this location from our forward model */ + { + double dd, refXYZ[3], delXYZ[3]; + fwddev(x, refXYZ, p->s[i].rgb); + if (vv < 0.5) + dd = 0.02; + else + dd = -0.02; + /* Matrix organization is J[XYZ][RGB] for del RGB->del XYZ*/ + for (j = 0; j < 3; j++) { + p->s[i].rgb[j] += dd; + fwddev(x, delXYZ, p->s[i].rgb); + p->s[i].j[0][j] = (delXYZ[0] - refXYZ[0]) / dd; + p->s[i].j[1][j] = (delXYZ[1] - refXYZ[1]) / dd; + p->s[i].j[2][j] = (delXYZ[2] - refXYZ[2]) / dd; + p->s[i].rgb[j] -= dd; + } + } +#endif + if (icmInverse3x3(p->s[i].ij, p->s[i].j)) { + error("dispcal: inverting Jacobian failed (2)"); + } + /* Make a copy of this Jacobian in case we get an invert failure later */ + icmCpy3x3(p->s[i].fb_ij, p->s[i].ij); + + /* Compute expected delta XYZ using new Jacobian */ + icmMulBy3x3(p->s[i].pdXYZ, p->s[i].j, p->s[i].pdrgb); + + p->s[i].de = b * os[j+1].de + (1.0 - b) * os[j].de; + p->s[i].dc = b * os[j+1].dc + (1.0 - b) * os[j].dc; + } + + free(os); +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +#ifdef NEVER +/* Do a linear interp of the ramdac */ +static void interp_ramdac(double cal[CAL_RES][3], double drgb[3], double rgb[3]) { + int i, j; + int gres = CAL_RES; + double w; + + /* For r,g & b */ + for (j = 0; j < 3; j++) { + int mi, gres_1 = gres-1; + double t, vv = rgb[j]; + t = gres * vv; + mi = (int)floor(t); /* Grid coordinate */ + if (mi < 0) /* Limit to valid cube base index range */ + mi = 0; + else if (mi >= gres_1) + mi = gres_1-1; + w = t - (double)mi; /* 1.0 - weight */ + + drgb[j] = (1.0 - w) * cal[mi][j] + w * cal[mi+1][j]; + } +} +#endif /* NEVER */ + +/* Given an XYZ, compute the color temperature and the delta E 2K to the locus */ +static double comp_ct( + double *de, /* If non-NULL, return CIEDE2000 to locus */ + double lxyz[3], /* If non-NULL, return normalised XYZ on locus */ + int plank, /* NZ if Plankian locus, 0 if Daylight locus */ + int dovct, /* NZ if visual match, 0 if traditional correlation */ + icxObserverType obType, /* If not default, set a custom observer */ + double xyz[3] /* Color to match */ +) { + double ct_xyz[3]; /* XYZ on locus */ + double nxyz[3]; /* Normalised input color */ + double ct, ctde; /* Color temperature & delta E to Black Body locus */ + icmXYZNumber wN; + + + if (obType == icxOT_default) + obType = icxOT_CIE_1931_2; + + if ((ct = icx_XYZ2ill_ct(ct_xyz, plank != 0 ? icxIT_Ptemp : icxIT_Dtemp, + obType, NULL, xyz, NULL, dovct)) < 0) + error ("Got bad color temperature conversion\n"); + + if (de != NULL) { + icmAry2XYZ(wN, ct_xyz); + icmAry2Ary(nxyz, xyz); + nxyz[0] /= xyz[1]; + nxyz[2] /= xyz[1]; + nxyz[1] /= xyz[1]; + ctde = icmXYZCIE2K(&wN, nxyz, ct_xyz); + *de = ctde; + } + if (lxyz != NULL) { + icmAry2Ary(lxyz, ct_xyz); + } + return ct; +} + +/* =================================================================== */ + +/* Default gamma */ +double g_def_gamma = 2.4; + +/* + + Flags used: + + ABCDEFGHIJKLMNOPQRSTUVWXYZ + upper .......... ....... . ... + lower ....... . ...... .... . + +*/ + +void usage(char *diag, ...) { + int i; + disppath **dp; + icompaths *icmps; + inst2_capability cap2 = inst2_none; + + fprintf(stderr,"Calibrate a Display, Version %s\n",ARGYLL_VERSION_STR); + fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n"); + if (setup_spyd2() == 2) + fprintf(stderr,"WARNING: This file contains a proprietary firmware image, and may not be freely distributed !\n"); + if (diag != NULL) { + va_list args; + fprintf(stderr,"Diagnostic: "); + va_start(args, diag); + vfprintf(stderr, diag, args); + va_end(args); + fprintf(stderr,"\n"); + } + fprintf(stderr,"usage: dispcal [options] outfile\n"); + fprintf(stderr," -v [n] Verbose mode\n"); +#if defined(UNIX_X11) + fprintf(stderr," -display displayname Choose X11 display name\n"); + fprintf(stderr," -d n[,m] Choose the display n from the following list (default 1)\n"); + fprintf(stderr," Optionally choose different display m for VideoLUT access\n"); +#else + fprintf(stderr," -d n Choose the display from the following list (default 1)\n"); +#endif + dp = get_displays(); + if (dp == NULL || dp[0] == NULL) + fprintf(stderr," ** No displays found **\n"); + else { + int i; + for (i = 0; ; i++) { + if (dp[i] == NULL) + break; + fprintf(stderr," %d = '%s'\n",i+1,dp[i]->description); + } + } + free_disppaths(dp); + fprintf(stderr," -dweb[:port] Display via a web server at port (default 8080)\n"); +// fprintf(stderr," -d fake Use a fake display device for testing, fake%s if present\n",ICC_FILE_EXT); + fprintf(stderr," -c listno Set communication port from the following list (default %d)\n",COMPORT); + if ((icmps = new_icompaths(g_log)) != NULL) { + icompath **paths; + if ((paths = icmps->paths) != NULL) { + int i; + for (i = 0; ; i++) { + if (paths[i] == NULL) + break; + if (paths[i]->itype == instSpyder2 && setup_spyd2() == 0) + fprintf(stderr," %d = '%s' !! Disabled - no firmware !!\n",i+1,paths[i]->name); + else + fprintf(stderr," %d = '%s'\n",i+1,paths[i]->name); + } + } else + fprintf(stderr," ** No ports found **\n"); + } + fprintf(stderr," -r Report on the calibrated display then exit\n"); + fprintf(stderr," -R Report on the uncalibrated display then exit\n"); + fprintf(stderr," -m Skip adjustment of the monitor controls\n"); + fprintf(stderr," -o [profile%s] Create fast matrix/shaper profile [different filename to outfile%s]\n",ICC_FILE_EXT,ICC_FILE_EXT); + fprintf(stderr," -O \"description\" Fast ICC Profile Description string (Default \"outfile\")\n"); + fprintf(stderr," -u Update previous calibration and (if -o used) ICC profile VideoLUTs\n"); + fprintf(stderr," -q [vlmh] Quality - Very Low, Low, Medium (def), High\n"); +// fprintf(stderr," -q [vfmsu] Speed - Very Fast, Fast, Medium (def), Slow, Ultra Slow\n"); + fprintf(stderr," -p Use telephoto mode (ie. for a projector) (if available)\n"); + cap2 = inst_show_disptype_options(stderr, " -y ", icmps, 0); + fprintf(stderr," -t [temp] White Daylight locus target, optional target temperaturee in deg. K (deflt.)\n"); + fprintf(stderr," -T [temp] White Black Body locus target, optional target temperaturee in deg. K\n"); + fprintf(stderr," -w x,y Set the target white point as chromaticity coordinates\n"); +#ifdef NEVER /* Not worth confusing people about this ? */ + fprintf(stderr," -L Show CCT/CDT rather than VCT/VDT during native white point adjustment\n"); +#endif + fprintf(stderr," -b bright Set the target white brightness in cd/m^2\n"); + fprintf(stderr," -g gamma Set the target response curve advertised gamma (Def. %3.1f)\n",g_def_gamma); + fprintf(stderr," Use \"-gl\" for L*a*b* curve\n"); + fprintf(stderr," Use \"-gs\" for sRGB curve\n"); + fprintf(stderr," Use \"-g709\" for REC 709 curve (should use -a as well!)\n"); + fprintf(stderr," Use \"-g240\" for SMPTE 240M curve (should use -a as well!)\n"); + fprintf(stderr," Use \"-G2.4 -f0\" for BT.1886\n"); + fprintf(stderr," -G gamma Set the target response curve actual technical gamma\n"); + fprintf(stderr," -f [degree] Amount of black level accounted for with output offset (default all output offset)\n"); + fprintf(stderr," -a ambient Use viewing condition adjustment for ambient in Lux\n"); + fprintf(stderr," -k factor Amount to correct black hue, 0 = none, 1 = full, Default = Automatic\n"); + fprintf(stderr," -A rate Rate of blending from neutral to black point. Default %.1f\n",NEUTRAL_BLEND_RATE); + fprintf(stderr," -B blkbright Set the target black brightness in cd/m^2\n"); + fprintf(stderr," -e [n] Run n verify passes on final curves\n"); + fprintf(stderr," -E Run only verify pass on installed calibration curves\n"); + fprintf(stderr," -P ho,vo,ss[,vs] Position test window and scale it\n"); + fprintf(stderr," ho,vi: 0.0 = left/top, 0.5 = center, 1.0 = right/bottom etc.\n"); + fprintf(stderr," ss: 0.5 = half, 1.0 = normal, 2.0 = double etc.\n"); + fprintf(stderr," -F Fill whole screen with black background\n"); +#if defined(UNIX_X11) + fprintf(stderr," -n Don't set override redirect on test window\n"); +#endif + fprintf(stderr," -J Run instrument calibration first (used rarely)\n"); + fprintf(stderr," -N Disable initial calibration of instrument if possible\n"); + fprintf(stderr," -H Use high resolution spectrum mode (if available)\n"); +// fprintf(stderr," -V Use adaptive measurement mode (if available)\n"); + if (cap2 & inst2_ccmx) + fprintf(stderr," -X file.ccmx Apply Colorimeter Correction Matrix\n"); + if (cap2 & inst2_ccss) { + fprintf(stderr," -X file.ccss Use Colorimeter Calibration Spectral Samples for calibration\n"); + fprintf(stderr," -Q observ Choose CIE Observer for spectrometer or CCSS colorimeter data:\n"); + fprintf(stderr," 1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2, 1964_10c\n"); + } + fprintf(stderr," -I b|w Drift compensation, Black: -Ib, White: -Iw, Both: -Ibw\n"); + fprintf(stderr," -Y A Use non-adaptive integration time mode (if available).\n"); + fprintf(stderr," -C \"command\" Invoke shell \"command\" each time a color is set\n"); + fprintf(stderr," -M \"command\" Invoke shell \"command\" each time a color is measured\n"); + fprintf(stderr," -W n|h|x Override serial port flow control: n = none, h = HW, x = Xon/Xoff\n"); + fprintf(stderr," -D [level] Print debug diagnostics to stderr\n"); + fprintf(stderr," inoutfile Base name for created or updated .cal and %s output files\n",ICC_FILE_EXT); + if (icmps != NULL) + icmps->del(icmps); + exit(1); +} + +int main(int argc, char *argv[]) { + int i, j, k; + int fa, nfa, mfa; /* current argument we're looking at */ + disppath *disp = NULL; /* Display being used */ + double hpatscale = 1.0, vpatscale = 1.0; /* scale factor for test patch size */ + double ho = 0.0, vo = 0.0; /* Test window offsets, -1.0 to 1.0 */ + int blackbg = 0; /* NZ if whole screen should be filled with black */ + int verb = 0; + int debug = 0; + int fake = 0; /* Use the fake device for testing */ + int override = 1; /* Override redirect on X11 */ + int docalib = 0; /* Do a manual instrument calibration */ + int doreport = 0; /* 1 = Report the current uncalibrated display response */ + /* 2 = Report the current calibrated display response */ + int docontrols = 1; /* Do adjustment of the display controls */ + int doprofile = 0; /* Create/update ICC profile */ + char *profDesc = NULL; /* Created profile description string */ + char *copyright = NULL; /* Copyright string */ + char *deviceMfgDesc = NULL; /* Device manufacturer string */ + char *modelDesc = NULL; /* Device model description string */ + int doupdate = 0; /* Do an update rather than a fresh calbration */ + int comport = COMPORT; /* COM port used */ + icompaths *icmps = NULL; + icompath *ipath = NULL; + flow_control fc = fc_nc; /* Default flow control */ + int dtype = 0; /* Display type selection charater */ + int tele = 0; /* nz if telephoto mode */ + int nocal = 0; /* Disable auto calibration */ + int highres = 0; /* Use high res mode if available */ + int nadaptive = 0; /* Use non-adaptive mode if available */ + int bdrift = 0; /* Flag, nz for black drift compensation */ + int wdrift = 0; /* Flag, nz for white drift compensation */ + double temp = 0.0; /* Color temperature (0 = native) */ + int planckian = 0; /* 0 = Daylight, 1 = Planckian color locus */ + int dovct = 1; /* Show VXT rather than CXT for adjusting white point */ + double wpx = 0.0, wpy = 0.0; /* White point xy (native) */ + double tbright = 0.0; /* Target white brightness ( 0.0 == max) */ + double gamma = 0.0; /* Advertised Gamma target */ + double egamma = 0.0; /* Effective Gamma target, NZ if set */ + double ambient = 0.0; /* NZ if viewing cond. adjustment to be used (cd/m^2) */ + double bkcorrect = -1.0; /* Level of black point correction, < 0 = auto */ + double bkbright = 0.0; /* Target black brightness ( 0.0 == min) */ + int quality = -99; /* Quality level, -2 = v, -1 = l, 0 = m, 1 = h, 2 = u */ + int isteps = 22; /* Initial measurement steps/3 (medium) */ + int rsteps = 64; /* Refinement measurement steps (medium) */ + double errthr = 1.5; /* Error threshold for refinement steps (medium) */ + int thrfail = 0; /* Set to NZ if failed to meet threshold target */ + double failerr = 0.0; /* Delta E of worst failed target */ + int mxits = 3; /* maximum iterations (medium) */ + int verify = 0; /* Do a verify after last refinement, 2 = do only verify. */ + int nver = 0; /* Number of verify passes after refinement */ + int webdisp = 0; /* NZ for web display, == port number */ + char *ccallout = NULL; /* Change color Shell callout */ + char *mcallout = NULL; /* Measure color Shell callout */ + char outname[MAXNAMEL+1] = { 0 }; /* Output cgats file base name */ + char iccoutname[MAXNAMEL+1] = { 0 };/* Output icc file base name */ + char ccxxname[MAXNAMEL+1] = "\000"; /* CCMX or CCSS file name */ + ccmx *cmx = NULL; /* Colorimeter Correction Matrix */ + ccss *ccs = NULL; /* Colorimeter Calibration Spectral Samples */ + int spec = 0; /* Want spectral data from instrument */ + icxObserverType obType = icxOT_default; + disprd *dr = NULL; /* Display patch read object */ + csamp asgrey; /* Main calibration loop test points */ + double dispLum = 0.0; /* Display luminence reading */ + int it; /* verify & refine iteration */ + int rv; + int fitord = 30; /* More seems to make curves smoother */ + int native = 1; /* 0 = use current or given calibration curve */ + /* 1 = set native linear op and use ramdac high prec'n */ + int noramdac = 0; /* Will be set to nz if can't set ramdac */ + int errc; /* Return value from new_disprd() */ + calx x; /* Context for calibration solution */ + + set_exe_path(argv[0]); /* Set global exe_path and error_program */ + check_if_not_interactive(); + +#if defined(__APPLE__) + { + SInt32 MacVers; + + /* Hmm. Maybe this should actually be 1.72 ?? */ + g_def_gamma = 1.8; + + /* OS X 10.6 uses a nominal gamma of 2.2 */ + if (Gestalt(gestaltSystemVersion, &MacVers) == noErr) { + if (MacVers >= 0x1060) { + g_def_gamma = 2.4; + } + } + } +#else + g_def_gamma = 2.4; /* Typical CRT gamma */ +#endif + gamma = g_def_gamma; + + setup_spyd2(); /* Load firware if available */ + + x.gammat = gt_power ; /* Default gamma type */ + x.egamma = 0.0; /* Default effective gamma none */ + x.oofff = 1.0; /* Default is all output ofset */ + x.vc = 0; /* No viewing conditions adjustment */ + x.svc = NULL; + x.dvc = NULL; + x.nbrate = NEUTRAL_BLEND_RATE; /* Rate of blending from black point to neutral axis */ + +#ifdef DEBUG_OFFSET + ho = 0.8; + vo = -0.8; +#endif + +#if defined(DEBUG) || defined(DEBUG_OFFSET) || defined(DEBUG_PLOT) + printf("!!!!!! Debug turned on !!!!!!\n"); +#endif + + if (argc <= 1) + usage("Too few arguments"); + + /* Process the arguments */ + mfa = 1; /* Minimum final arguments */ + for(fa = 1;fa < argc;fa++) { + nfa = fa; /* skip to nfa if next argument is used */ + if (argv[fa][0] == '-') /* Look for any flags */ + { + char *na = NULL; /* next argument after flag, null if none */ + + if (argv[fa][2] != '\000') + na = &argv[fa][2]; /* next is directly after flag */ + else { + if ((fa+1+mfa) < argc) { + if (argv[fa+1][0] != '-') + { + nfa = fa + 1; + na = argv[nfa]; /* next is seperate non-flag argument */ + } + } + } + + if (argv[fa][1] == '?' || argv[fa][1] == '-') { + usage("Usage requested"); + + } else if (argv[fa][1] == 'v') { + verb = 1; + if (na != NULL && na[0] >= '0' && na[0] <= '9') { + verb = atoi(na); + fa = nfa; + } + g_log->verb = verb; + + /* Display number */ + } else if (argv[fa][1] == 'd') { + if (strncmp(na,"web",3) == 0 + || strncmp(na,"WEB",3) == 0) { + webdisp = 8080; + if (na[3] == ':') { + webdisp = atoi(na+4); + if (webdisp == 0 || webdisp > 65535) + usage("Web port number must be in range 1..65535"); + } + fa = nfa; + } else { +#if defined(UNIX_X11) + int ix, iv; + + if (strcmp(&argv[fa][2], "isplay") == 0 || strcmp(&argv[fa][2], "ISPLAY") == 0) { + if (++fa >= argc || argv[fa][0] == '-') usage("Parameter expected following -display"); + setenv("DISPLAY", argv[fa], 1); + } else { + if (na == NULL) usage("Parameter expected following -d"); + fa = nfa; + if (strcmp(na,"fake") == 0) { + fake = 1; + } else { + if (sscanf(na, "%d,%d",&ix,&iv) != 2) { + ix = atoi(na); + iv = 0; + } + if (disp != NULL) + free_a_disppath(disp); + if ((disp = get_a_display(ix-1)) == NULL) + usage("-d parameter %d out of range",ix); + if (iv > 0) + disp->rscreen = iv-1; + } + } +#else + int ix; + if (na == NULL) usage("Parameter expected following -d"); + fa = nfa; + if (strcmp(na,"fake") == 0) { + fake = 1; + } else { + ix = atoi(na); + if (disp != NULL) + free_a_disppath(disp); + if ((disp = get_a_display(ix-1)) == NULL) + usage("-d parameter %d out of range",ix); + } +#endif + } + + } else if (argv[fa][1] == 'J') { + docalib = 1; + + } else if (argv[fa][1] == 'N') { + nocal = 1; + + /* High res mode */ + } else if (argv[fa][1] == 'H') { + highres = 1; + + /* Adaptive mode - now default, so flag is deprecated */ + } else if (argv[fa][1] == 'V') { + warning("dispcal -V flag is deprecated"); + + /* Colorimeter Correction Matrix */ + /* or Colorimeter Calibration Spectral Samples */ + } else if (argv[fa][1] == 'X') { + int ix; + fa = nfa; + if (na == NULL) usage("Parameter expected following -X"); + strncpy(ccxxname,na,MAXNAMEL-1); ccxxname[MAXNAMEL-1] = '\000'; + + /* Drift Compensation */ + } else if (argv[fa][1] == 'I') { + fa = nfa; + if (na == NULL || na[0] == '\000') usage("Parameter expected after -I"); + for (i=0; ; i++) { + if (na[i] == '\000') + break; + if (na[i] == 'b' || na[i] == 'B') + bdrift = 1; + else if (na[i] == 'w' || na[i] == 'W') + wdrift = 1; + else + usage("-I parameter '%c' not recognised",na[i]); + } + + /* Spectral Observer type */ + } else if (argv[fa][1] == 'Q') { + fa = nfa; + if (na == NULL) usage("Parameter expecte after -Q"); + if (strcmp(na, "1931_2") == 0) { /* Classic 2 degree */ + obType = icxOT_CIE_1931_2; + } else if (strcmp(na, "1964_10") == 0) { /* Classic 10 degree */ + obType = icxOT_CIE_1964_10; + } else if (strcmp(na, "1964_10c") == 0) { /* 10 degree corrected */ + obType = icxOT_CIE_1964_10c; + } else if (strcmp(na, "1955_2") == 0) { /* Stiles and Burch 1955 2 degree */ + obType = icxOT_Stiles_Burch_2; + } else if (strcmp(na, "1978_2") == 0) { /* Judd and Voss 1978 2 degree */ + obType = icxOT_Judd_Voss_2; + } else if (strcmp(na, "shaw") == 0) { /* Shaw and Fairchilds 1997 2 degree */ + obType = icxOT_Shaw_Fairchild_2; + } else + usage("-Q parameter '%s' not recognised",na); + + /* Change color callout */ + } else if (argv[fa][1] == 'C') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -C"); + ccallout = na; + + /* Measure color callout */ + } else if (argv[fa][1] == 'M') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -M"); + mcallout = na; + + /* Serial port flow control */ + } else if (argv[fa][1] == 'W') { + fa = nfa; + if (na == NULL) usage("Paramater expected following -W"); + if (na[0] == 'n' || na[0] == 'N') + fc = fc_none; + else if (na[0] == 'h' || na[0] == 'H') + fc = fc_Hardware; + else if (na[0] == 'x' || na[0] == 'X') + fc = fc_XonXOff; + else + usage("-W parameter '%c' not recognised",na[0]); + + /* Debug coms */ + } else if (argv[fa][1] == 'D') { + debug = 1; + if (na != NULL && na[0] >= '0' && na[0] <= '9') { + debug = atoi(na); + fa = nfa; + } + g_log->debug = debug; + callback_ddebug = 1; /* dispwin global */ + + /* Black point correction amount */ + } else if (argv[fa][1] == 'k') { + fa = nfa; + if (na == NULL) usage("Paramater expected following -k"); + bkcorrect = atof(na); + if (bkcorrect < 0.0 || bkcorrect > 1.0) usage("-k parameter must be between 0.0 and 1.0"); + /* Neutral blend rate (power) */ + } else if (argv[fa][1] == 'A') { + fa = nfa; + if (na == NULL) usage("Paramater expected following -A"); + x.nbrate = atof(na); + if (x.nbrate < 0.05 || x.nbrate > 20.0) usage("-A parameter must be between 0.05 and 20.0"); + /* Black brightness */ + } else if (argv[fa][1] == 'B') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -B"); + bkbright = atof(na); + if (bkbright <= 0.0 || bkbright > 100000.0) usage("-B parameter %f out of range",bkbright); + + /* Number of verify passes */ + } else if (argv[fa][1] == 'e') { + verify = 1; + nver = 1; + if (na != NULL && na[0] >= '0' && na[0] <= '9') { + nver = atoi(na); + fa = nfa; + } + + } else if (argv[fa][1] == 'E') { + verify = 2; + mfa = 0; + +#if defined(UNIX_X11) + } else if (argv[fa][1] == 'n') { + override = 0; +#endif /* UNIX */ + /* COM port */ + } else if (argv[fa][1] == 'c') { + fa = nfa; + if (na == NULL) usage("Paramater expected following -c"); + comport = atoi(na); + if (comport < 1 || comport > 50) usage("-c parameter %d out of range",comport); + + /* Telephoto */ + } else if (argv[fa][1] == 'p') { + tele = 1; + + } else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') { + if (argv[fa][1] == 'R') + doreport = 1; /* raw */ + else + doreport = 2; /* Calibrated */ + mfa = 0; + + } else if (argv[fa][1] == 'm') { + docontrols = 0; + + /* Output/update ICC profile [optional different name] */ + } else if (argv[fa][1] == 'o') { + doprofile = 1; + + if (na != NULL) { /* Found an optional icc profile name */ + fa = nfa; + strncpy(iccoutname,na,MAXNAMEL); iccoutname[MAXNAMEL] = '\000'; + } + + /* Fast Profile Description */ + } else if (argv[fa][1] == 'O') { + fa = nfa; + if (na == NULL) usage("Expect argument to profile description flag -O"); + profDesc = na; + + /* Update calibration and (optionally) profile */ + } else if (argv[fa][1] == 'u') { + doupdate = 1; + docontrols = 0; + + /* Speed/Quality */ + } else if (argv[fa][1] == 'q') { + fa = nfa; + if (na == NULL) usage("Parameter expected following -q"); + switch (na[0]) { + case 'L': /* Test value */ + quality = -3; + break; + case 'v': /* very fast */ + quality = -2; + break; + case 'f': /* fast */ + case 'l': + quality = -1; + break; + case 'm': /* medium */ + case 'M': + quality = 0; + break; + case 's': /* slow */ + case 'h': + case 'H': + quality = 1; + break; + case 'u': /* ultra slow */ + case 'U': + quality = 2; + break; + default: + usage("-q parameter '%c' not recognised",na[0]); + } + + /* Display type */ + } else if (argv[fa][1] == 'y') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -y"); + dtype = na[0]; + + /* Daylight color temperature */ + } else if (argv[fa][1] == 't' || argv[fa][1] == 'T') { + if (argv[fa][1] == 'T') + planckian = 1; + else + planckian = 0; + if (na != NULL) { + fa = nfa; + temp = atof(na); + if (temp < 1000.0 || temp > 15000.0) usage("-%c parameter %f out of range",argv[fa][1], temp); + } + + /* White point as x, y */ + } else if (argv[fa][1] == 'w') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -w"); + if (sscanf(na, " %lf,%lf ", &wpx, &wpy) != 2) + usage("-w parameter '%s' not recognised",na); + + /* Show CXT rather than VXT when adjusting native white point */ + } else if (argv[fa][1] == 'L') { + dovct = 0; + + /* White brightness */ + } else if (argv[fa][1] == 'b') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -b"); + tbright = atof(na); + if (tbright <= 0.0 || tbright > 100000.0) usage("-b parameter %f out of range",tbright); + + /* Target transfer curve */ + } else if (argv[fa][1] == 'g') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -g"); + if ((na[0] == 'l' || na[0] == 'L') && na[1] == '\000') + x.gammat = gt_Lab; + else if ((na[0] == 's' || na[0] == 'S') && na[1] == '\000') + x.gammat = gt_sRGB; + else if (strcmp(na, "709") == 0) + x.gammat = gt_Rec709; + else if (strcmp(na, "240") == 0) + x.gammat = gt_SMPTE240M; + else { + gamma = atof(na); + if (gamma <= 0.0 || gamma > 10.0) usage("-g parameter %f out of range",gamma); + x.gammat = gt_power; + } + + /* Effective gamma power */ + } else if (argv[fa][1] == 'G') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -G"); + egamma = atof(na); + if (egamma <= 0.0 || egamma > 10.0) usage("-G parameter %f out of range",egamma); + x.gammat = gt_power; + + /* Degree of output offset */ + } else if (argv[fa][1] == 'f') { + fa = nfa; + if (na == NULL) { + x.oofff = 0.0; + } else { + x.oofff = atof(na); + if (x.oofff < 0.0 || x.oofff > 1.0) + usage("-f parameter %f out of range",x.oofff); + } + + /* Ambient light level */ + } else if (argv[fa][1] == 'a') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -a"); + ambient = atof(na); + if (ambient < 0.0) + usage("-a parameter %f out of range",ambient); + ambient /= 3.141592654; /* Convert from Lux to cd/m^2 */ + + /* Test patch offset and size */ + } else if (argv[fa][1] == 'P') { + fa = nfa; + if (na == NULL) usage("Parameter expected after -P"); + if (sscanf(na, " %lf,%lf,%lf,%lf ", &ho, &vo, &hpatscale, &vpatscale) == 4) { + ; + } else if (sscanf(na, " %lf,%lf,%lf ", &ho, &vo, &hpatscale) == 3) { + vpatscale = hpatscale; + } else { + usage("-P parameter '%s' not recognised",na); + } + if (ho < 0.0 || ho > 1.0 + || vo < 0.0 || vo > 1.0 + || hpatscale <= 0.0 || hpatscale > 50.0 + || vpatscale <= 0.0 || vpatscale > 50.0) + usage("-P parameters %f %f %f %f out of range",ho,vo,hpatscale,vpatscale); + ho = 2.0 * ho - 1.0; + vo = 2.0 * vo - 1.0; + + /* Black background */ + } else if (argv[fa][1] == 'F') { + blackbg = 1; + + /* Extra flags */ + } else if (argv[fa][1] == 'Y') { + if (na == NULL) + usage("Flag '-Y' expects extra flag"); + + if (na[0] == 'A') { + nadaptive = 1; + } else { + usage("Flag '-Y %c' not recognised",na[0]); + } + + } else + usage("Flag '-%c' not recognised",argv[fa][1]); + } else + break; + } + + /* No explicit display has been set */ + if ( +#ifndef SHOW_WINDOW_ONFAKE + !fake && +#endif + disp == NULL) { + int ix = 0; +#if defined(UNIX_X11) + char *dn, *pp; + + if ((dn = getenv("DISPLAY")) != NULL) { + if ((pp = strrchr(dn, ':')) != NULL) { + if ((pp = strchr(pp, '.')) != NULL) { + if (pp[1] != '\000') + ix = atoi(pp+1); + } + } + } +#endif + if ((disp = get_a_display(ix)) == NULL) + error("Unable to open the default display"); + } + + /* See if there is an environment variable ccxx */ + if (ccxxname[0] == '\000') { + char *na; + if ((na = getenv("ARGYLL_COLMTER_CAL_SPEC_SET")) != NULL) { + strncpy(ccxxname,na,MAXNAMEL-1); ccxxname[MAXNAMEL-1] = '\000'; + + } else if ((na = getenv("ARGYLL_COLMTER_COR_MATRIX")) != NULL) { + strncpy(ccxxname,na,MAXNAMEL-1); ccxxname[MAXNAMEL-1] = '\000'; + } + } + + /* Load up CCMX or CCSS */ + if (ccxxname[0] != '\000') { + if ((cmx = new_ccmx()) == NULL + || cmx->read_ccmx(cmx, ccxxname)) { + if (cmx != NULL) { + cmx->del(cmx); + cmx = NULL; + } + + /* CCMX failed, try CCSS */ + if ((ccs = new_ccss()) == NULL + || ccs->read_ccss(ccs, ccxxname)) { + if (ccs != NULL) { + ccs->del(ccs); + ccs = NULL; + error("Reading CCMX/CCSS File '%s' failed\n", ccxxname); + } + } + } + } + + if (fake) + comport = -99; + if ((icmps = new_icompaths(g_log)) == NULL) + error("Finding instrument paths failed"); + if ((ipath = icmps->get_path(icmps, comport)) == NULL) + error("No instrument at port %d",comport); + + if (docalib) { + if ((rv = disprd_calibration(ipath, fc, dtype, 0, tele, nadaptive, nocal, disp, + webdisp, blackbg, override, + 100.0 * hpatscale, 100.0 * vpatscale, + ho, vo, g_log)) != 0) { + error("docalibration failed with return value %d\n",rv); + } + } + + if (verify != 2 && doreport == 0) { + /* Get the file name argument */ + if (fa >= argc || argv[fa][0] == '-') usage("Output filname parameter not found"); + strncpy(outname,argv[fa],MAXNAMEL-4); outname[MAXNAMEL-4] = '\000'; + strcat(outname,".cal"); + if (iccoutname[0] == '\000') { + strncpy(iccoutname,argv[fa++],MAXNAMEL-4); iccoutname[MAXNAMEL-4] = '\000'; + strcat(iccoutname,ICC_FILE_EXT); + } + } + + if (verify == 2) { + if (doupdate) + warning("Update flag ignored because we're doing a verify only"); + doupdate = 0; + docontrols = 0; + } + + if (doreport != 0) { + if (verify == 2) + warning("Verify flag ignored because we're doing a report only"); + verify = 0; + } + + /* Normally calibrate against native response */ + if (verify == 2 || doreport == 2) + native = 0; /* But measure calibrated response of verify or report calibrated */ + + /* Get ready to do some readings */ + if ((dr = new_disprd(&errc, ipath, fc, dtype, 0, tele, nadaptive, nocal, + highres, native, &noramdac, NULL, 0, 0, disp, blackbg, override, + webdisp, ccallout, mcallout, + 100.0 * hpatscale, 100.0 * vpatscale, ho, vo, + cmx != NULL ? cmx->matrix : NULL, + ccs != NULL ? ccs->samples : NULL, ccs != NULL ? ccs->no_samp : 0, + spec, obType, NULL, bdrift, wdrift, + "fake" ICC_FILE_EXT, g_log)) == NULL) + error("new_disprd() failed with '%s'\n",disprd_err(errc)); + + if (native != 0 && noramdac) { + warning("No access to VideoLUTs, so calibrating display as-is rather than native"); + if (doprofile) + warning("Profile will reflect the as-is display response and not contain a 'vcgt' tag"); + native = 0; + + if (doupdate && doprofile) + error("Can't update a profile that doesn't use the 'vcgt' tag for calibration"); + } + + if (icmps != NULL) { + icmps->del(icmps); + icmps = NULL; + } + if (cmx != NULL) { + cmx->del(cmx); + cmx = NULL; + } + if (ccs != NULL) { + ccs->del(ccs); + ccs = NULL; + } + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + if (doreport) { + col tcols[3] = { /* Base set of test colors */ + { 0.0, 0.0, 0.0 }, + { 0.5, 0.5, 0.5 }, + { 1.0, 1.0, 1.0 } + }; + double cct, cct_de; /* Color temperatures and DE 2K */ + double cdt, cdt_de; + double vct, vct_de; + double vdt, vdt_de; + double cgamma, w[3], wp[2]; + int sigbits = 0; /* Number of significant bits in VideoLUT/display/instrument */ + + if ((rv = dr->read(dr, tcols, 3, 1, 3, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + +//printf("~1 Got black = %f, half = %f, white = %f\n",tcols[0].XYZ[1],tcols[1].XYZ[1],tcols[2].XYZ[1]); + /* Normalised XYZ white point */ + w[0] = tcols[2].XYZ[0]/tcols[2].XYZ[1]; + w[1] = tcols[2].XYZ[1]/tcols[2].XYZ[1]; + w[2] = tcols[2].XYZ[2]/tcols[2].XYZ[1]; + + /* White point chromaticity coordinates */ + wp[0] = w[0]/(w[0] + w[1] + w[2]); + wp[1] = w[1]/(w[0] + w[1] + w[2]); + + cct = comp_ct(&cct_de, NULL, 1, 0, obType, w); /* Compute CCT */ + cdt = comp_ct(&cdt_de, NULL, 0, 0, obType, w); /* Compute CDT */ + vct = comp_ct(&vct_de, NULL, 1, 1, obType, w); /* Compute VCT */ + vdt = comp_ct(&vdt_de, NULL, 0, 1, obType, w); /* Compute VDT */ + + /* Compute advertised current gamma - use the gross curve shape for robustness */ + cgamma = pop_gamma(tcols[0].XYZ[1], tcols[1].XYZ[1], tcols[2].XYZ[1]); + +#ifdef MEAS_RES + /* See if we can detect what sort of precision the LUT entries */ + /* have. Our ability to detect this may be limited by the instrument */ + /* (ie. Huey and Spyder 2) */ + if (doreport == 1) { +#define MAX_RES_SAMPS 24 + col ttt[MAX_RES_SAMPS]; + int res_samps = 6; + double a0, a1, a2, dd; + int n; + int issig = 0; + + if (verb) + printf("Measuring VideoLUT table entry precision.\n"); + + /* Run a small state machine until we come to a conclusion */ + sigbits = 8; + for (issig = 0; issig < 2; ) { + + DBG((dbgo,"Trying %d bits\n",sigbits)); + /* Do the test */ + for (n = 0; n < res_samps; n++) { + double v; +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(sigbits); +#endif + /* Notional test value */ + v = (5 << (sigbits-3))/((1 << sigbits) - 1.0); + /* And -1, 0 , +1 bit test values */ + if ((n % 3) == 2) + v += 1.0/((1 << sigbits) - 1.0); + else if ((n % 3) == 1) + v += 0.0/((1 << sigbits) - 1.0); + else + v += -1.0/((1 << sigbits) - 1.0); + ttt[n].r = ttt[n].g = ttt[n].b = v; + } + if ((rv = dr->read(dr, ttt, res_samps, 1, res_samps, 1, 0, instNoClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + /* Average the readings for each test value */ + a0 = a1 = a2 = 0.0; + for (n = 0; n < res_samps; n++) { + double v = ttt[n].XYZ[1]; + if ((n % 3) == 2) { + a2 += v; + } else if ((n % 3) == 1) { + a1 += v; + } else { + a0 += v; + } + } + a0 /= (res_samps / 3.0); + a1 /= (res_samps / 3.0); + a2 /= (res_samps / 3.0); + /* Judge significance of any differences */ + dd = 0.0; + for (n = 0; n < res_samps; n++) { + double tt; + if ((n % 3) == 2) + tt = fabs(a2 - ttt[n].XYZ[1]); + else if ((n % 3) == 1) + tt = fabs(a1 - ttt[n].XYZ[1]); + else + tt = fabs(a0 - ttt[n].XYZ[1]); + dd += tt * tt; + } + dd /= res_samps; + dd = sqrt(dd); + if (fabs(a1 - a0) > (2.0 * dd) && fabs(a2 - a1) > (2.0 * dd)) + issig = 1; /* Noticable difference */ + else + issig = 0; /* No noticable difference */ + DBG((dbgo,"Bits %d: Between = %f, %f within = %f, sig = %s\n",sigbits, fabs(a1 - a0), fabs(a2 - a1),dd, issig ? "yes" : "no")); + + switch(sigbits) { + case 8: /* Do another trial */ + if (issig) { + sigbits = 10; + res_samps = 6; + } else { + sigbits = 6; + } + break; + case 6: /* Do another trial or give up */ + if (issig) { + sigbits = 7; + } else { + sigbits = 0; + issig = 2; /* Give up */ + } + break; + case 7: /* Terminal */ + if (!issig) + sigbits = 6; + issig = 2; /* Stop here */ + break; + case 10: /* Do another trial */ + if (issig) { + sigbits = 12; + res_samps = 9; + } else { + sigbits = 9; + } + break; + case 12: /* Do another trial or give up */ + if (issig) { + issig = 2; /* Stop here */ + } else { + sigbits = 11; + } + break; + case 11: /* Terminal */ + if (!issig) + sigbits = 10; + issig = 2; /* Stop here */ + break; + case 9: /* Terminal */ + if (!issig) + sigbits = 8; + issig = 2; /* Stop here */ + break; + + default: + error("Unexpected number of bits in depth test (bits)",sigbits); + } + } + } +#endif /* MEAS_RES */ + + if (doreport == 2) + printf("Current calibration response:\n"); + else + printf("Uncalibrated response:\n"); + printf("Black level = %.2f cd/m^2\n",tcols[0].XYZ[1]); + printf("White level = %.2f cd/m^2\n",tcols[2].XYZ[1]); + printf("Aprox. gamma = %.2f\n",cgamma); + printf("Contrast ratio = %.0f:1\n",tcols[2].XYZ[1]/tcols[0].XYZ[1]); + printf("White chromaticity coordinates %.4f, %.4f\n",wp[0],wp[1]); + printf("White Correlated Color Temperature = %.0fK, DE 2K to locus = %4.1f\n",cct,cct_de); + printf("White Correlated Daylight Temperature = %.0fK, DE 2K to locus = %4.1f\n",cdt,cdt_de); + printf("White Visual Color Temperature = %.0fK, DE 2K to locus = %4.1f\n",vct,vct_de); + printf("White Visual Daylight Temperature = %.0fK, DE 2K to locus = %4.1f\n",vdt,vdt_de); +#ifdef MEAS_RES + if (doreport == 1) { + if (sigbits == 0) { + warning("Unable to determine LUT entry bit depth - problems ?"); + } else { + printf("Effective LUT entry depth seems to be %d bits\n",sigbits); + } + } +#endif /* MEAS_RES */ + dr->del(dr); + exit(0); + } + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* If we're updating, retrieve the previously used settings, */ + /* and device model */ + if (doupdate) { + cgats *icg; /* output cgats structure */ + int nsamp; + mcvco *rdv[3]; /* Scattered data for ramdac curves */ + int fi; /* Field index */ + int si[4]; /* Set fields */ + + if (verb) { + if (doprofile) + printf("Updating previous calibration and profile\n"); + else + printf("Updating previous calibration\n"); + } + + icg = new_cgats(); /* Create a CGATS structure */ + icg->add_other(icg, "CAL"); /* our special type is Calibration file */ + + if (icg->read_name(icg, outname)) { + dr->del(dr); + error("Can't update '%s' - read error : %s",outname, icg->err); + } + + if (icg->ntables == 0 + || icg->t[0].tt != tt_other || icg->t[0].oi != 0 + || icg->t[1].tt != tt_other || icg->t[1].oi != 0) { + dr->del(dr); + error("Can't update '%s' - wrong type of file",outname); + } + if (icg->ntables < 2) { + dr->del(dr); + error("Can't update '%s' - there aren't two tables",outname); + } + +//printf("~1 reading previous cal, got 2 tables\n"); + + /* Read in the setup, user and model values */ + + if (dtype == 0) { /* If the use hasn't set anything */ + if ((fi = icg->find_kword(icg, 0, "DEVICE_TYPE")) >= 0) { + if (strcmp(icg->t[0].kdata[fi], "CRT") == 0) + dtype = 'c'; + else if (strcmp(icg->t[0].kdata[fi], "LCD") == 0) + dtype = 'l'; + else + dtype = icg->t[0].kdata[fi][0]; + } + } +//printf("~1 dealt with device type\n"); + + if ((fi = icg->find_kword(icg, 0, "TARGET_WHITE_XYZ")) < 0) { + dr->del(dr); + error ("Can't update '%s' - can't find field 'TARGET_WHITE_XYZ'",outname); + } + if (sscanf(icg->t[0].kdata[fi], "%lf %lf %lf", &x.twh[0], &x.twh[1], &x.twh[2]) != 3) { + dr->del(dr); + error ("Can't update '%s' - reading field 'TARGET_WHITE_XYZ' failed",outname); + } + x.nwh[0] = x.twh[0] / x.twh[1]; + x.nwh[1] = x.twh[1] / x.twh[1]; + x.nwh[2] = x.twh[2] / x.twh[1]; + + if ((fi = icg->find_kword(icg, 0, "NATIVE_TARGET_WHITE")) >= 0) { + wpx = wpy = 0.0; + temp = 0.0; + tbright = 0.0; + } else { + wpx = wpy = 0.0001; + temp = 1.0; + tbright = 1.0; + } +//printf("~1 dealt with target white\n"); + + if ((fi = icg->find_kword(icg, 0, "TARGET_GAMMA")) < 0) { + dr->del(dr); + error ("Can't update '%s' - can't find field 'TARGET_GAMMA'",outname); + } + if (strcmp(icg->t[0].kdata[fi], "L_STAR") == 0) + x.gammat = gt_Lab; + else if (strcmp(icg->t[0].kdata[fi], "sRGB") == 0) + x.gammat = gt_sRGB; + else if (strcmp(icg->t[0].kdata[fi], "REC709") == 0) + x.gammat = gt_Rec709; + else if (strcmp(icg->t[0].kdata[fi], "SMPTE240M") == 0) + x.gammat = gt_SMPTE240M; + else { + x.gammat = 0; + gamma = atof(icg->t[0].kdata[fi]); + + if (fabs(gamma) < 0.1 || fabs(gamma) > 5.0) { + dr->del(dr); + error ("Can't update '%s' - field 'TARGET_GAMMA' has bad value %f",outname,fabs(gamma)); + } + if (gamma < 0.0) { /* Effective gamma = actual power value */ + egamma = -gamma; + } + } + if ((fi = icg->find_kword(icg, 0, "DEGREE_OF_BLACK_OUTPUT_OFFSET")) < 0) { + /* Backward compatibility if value is not present */ + if (x.gammat == gt_Lab || x.gammat == gt_sRGB) + x.oofff = 1.0; + else + x.oofff = 0.0; + } else { + x.oofff = atof(icg->t[0].kdata[fi]); + } + if ((fi = icg->find_kword(icg, 0, "TARGET_BLACK_BRIGHTNESS")) < 0) { + bkbright = 0.0; /* Native */ + } else { + bkbright = atof(icg->t[0].kdata[fi]); + } + + + if ((fi = icg->find_kword(icg, 0, "BLACK_POINT_CORRECTION")) < 0) { + dr->del(dr); + error ("Can't update '%s' - can't find field 'BLACK_POINT_CORRECTION'",outname); + } + bkcorrect = atof(icg->t[0].kdata[fi]); + if (bkcorrect < 0.0 || bkcorrect > 1.0) { + dr->del(dr); + error ("Can't update '%s' - field 'BLACK_POINT_CORRECTION' has bad value %f",outname,bkcorrect); + } + + if ((fi = icg->find_kword(icg, 0, "BLACK_NEUTRAL_BLEND_RATE")) < 0) { + x.nbrate = 8.0; /* Backwards compatibility value */ + } else { + x.nbrate = atof(icg->t[0].kdata[fi]); + if (x.nbrate < 0.05 || x.nbrate > 20.0) { + dr->del(dr); + error ("Can't update '%s' - field 'BLACK_NEUTRAL_BLEND_RATE' has bad value %f",outname,x.nbrate); + } + } + + if ((fi = icg->find_kword(icg, 0, "QUALITY")) < 0) { + dr->del(dr); + error ("Can't update '%s' - can't find field 'QUALITY'",outname); + } + if (quality < -50) { /* User hasn't overridden quality */ + if (strcmp(icg->t[0].kdata[fi], "ultra low") == 0) + quality = -3; + else if (strcmp(icg->t[0].kdata[fi], "very low") == 0) + quality = -2; + else if (strcmp(icg->t[0].kdata[fi], "low") == 0) + quality = -1; + else if (strcmp(icg->t[0].kdata[fi], "medium") == 0) + quality = 0; + else if (strcmp(icg->t[0].kdata[fi], "high") == 0) + quality = 1; + else if (strcmp(icg->t[0].kdata[fi], "ultra high") == 0) + quality = 2; + else { + dr->del(dr); + error ("Can't update '%s' - field 'QUALITY' has unrecognised value '%s'", + outname,icg->t[0].kdata[fi]); + } + } +//printf("~1 dealt with quality\n"); + + /* Read in the last set of calibration curves used */ + if ((nsamp = icg->t[0].nsets) < 2) { + dr->del(dr); + error("Can't update '%s' - %d not enough data points in calibration curves", + outname,nsamp); + } +//printf("~1 got %d points in calibration curves\n",nsamp); + + for (k = 0; k < 3; k++) { + if ((x.rdac[k] = new_mcv()) == NULL) { + dr->del(dr); + error("new_mcv x.rdac[%d] failed",k); + } + if ((rdv[k] = malloc(sizeof(mcvco) * nsamp)) == NULL) { + dr->del(dr); + error ("Malloc of scattered data points failed"); + } + } +//printf("~1 allocated calibration curve objects\n"); + + /* Read the current calibration curve points (usually CAL_RES of them) */ + for (k = 0; k < 4; k++) { + char *fnames[4] = { "RGB_I", "RGB_R", "RGB_G", "RGB_B" }; + + if ((si[k] = icg->find_field(icg, 0, fnames[k])) < 0) { + dr->del(dr); + error ("Can't updata '%s' - can't find field '%s'",outname,fnames[k]); + } + if (icg->t[0].ftype[si[k]] != r_t) { + dr->del(dr); + error ("Can't updata '%s' - field '%s' is wrong type",outname,fnames[k]); + } + } +//printf("~1 Found calibration curve fields\n"); + + for (i = 0; i < nsamp; i++) { + rdv[0][i].p = + rdv[1][i].p = + rdv[2][i].p = + *((double *)icg->t[0].fdata[i][si[0]]); + for (k = 0; k < 3; k++) { /* RGB */ + rdv[k][i].v = *((double *)icg->t[0].fdata[i][si[k + 1]]); + } + rdv[0][i].w = rdv[1][i].w = rdv[2][i].w = 1.0; + } +//printf("~1 Read calibration curve data points\n"); + for (k = 0; k < 3; k++) { + x.rdac[k]->fit(x.rdac[k], 0, fitord, rdv[k], nsamp, RDAC_SMOOTH); + free (rdv[k]); + } +//printf("~1 Fitted calibration curves\n"); + + /* Read in the per channel forward model curves */ + for (k = 0; k < 3; k++) { + char *fnames[3] = { "R_P", "G_P", "B_P" }; + double *pp; + +//printf("~1 Reading device curve channel %d\n",k); + if ((si[k] = icg->find_field(icg, 1, fnames[k])) < 0) { + dr->del(dr); + error ("Can't updata '%s' - can't find field '%s'",outname,fnames[k]); + } + if (icg->t[1].ftype[si[k]] != r_t) { + dr->del(dr); + error ("Can't updata '%s' - field '%s' is wrong type",outname,fnames[k]); + } + /* Create the model curves */ + if ((pp = (double *)malloc(icg->t[1].nsets * sizeof(double))) == NULL) { + dr->del(dr); + error ("Malloc of device curve parameters"); + } + for (i = 0; i < icg->t[1].nsets; i++) + pp[i] = *((double *)icg->t[1].fdata[i][si[k]]); + + if ((x.dcvs[k] = new_mcv_p(pp, icg->t[1].nsets)) == NULL) { + dr->del(dr); + error("new_mcv x.dcvs[%d] failed",k); + } + free(pp); + } + + icg->del(icg); +//printf("~1 read in previous settings and device model\n"); + } + + /* Be nice - check we can read the iccprofile before calibrating the display */ + if (verify != 2 && doupdate && doprofile) { + icmFile *ic_fp; + icc *icco; + + if ((icco = new_icc()) == NULL) { + dr->del(dr); + error ("Creation of ICC object to read profile '%s' failed",iccoutname); + } + + /* Open up the profile for reading */ + if ((ic_fp = new_icmFileStd_name(iccoutname,"r")) == NULL) { + dr->del(dr); + error ("Can't open file '%s'",iccoutname); + } + + /* Read header etc. */ + if ((rv = icco->read(icco,ic_fp,0)) != 0) { + dr->del(dr); + error ("Reading profile '%s' failed with %d, %s",iccoutname, rv,icco->err); + } + + ic_fp->del(ic_fp); + + if (icco->find_tag(icco, icSigVideoCardGammaTag) != 0) { + dr->del(dr); + error("Can't find VideoCardGamma tag in file '%s': %d, %s", + iccoutname, icco->errc,icco->err); + } + icco->del(icco); + } + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* Convert quality level to iterations etc. */ + /* Note that final tolerance is often double the */ + /* final errth, because one more corrections is always */ + /* performed after the last reading. */ + switch (quality) { + case -3: /* Test value */ + isteps = 3; + rsteps = 9; + mxits = 1; + errthr = 2.0; + break; + case -2: /* Very low */ + isteps = 10; + rsteps = 16; + errthr = 1.5; + if (doupdate) + mxits = 1; + else + mxits = 1; + break; + case -1: /* Low */ + if (verify != 2 && doprofile && !doupdate) + isteps = 24; /* Use more steps if we're creating a profile */ + else + isteps = 12; + rsteps = 32; + errthr = 0.9; + if (doupdate) + mxits = 1; + else + mxits = 2; + break; + default: + case 0: /* Medum */ + quality = 0; /* In case it wasn't set */ + if (verify != 2 && doprofile && !doupdate) + isteps = 32; /* Use more steps if we're creating a profile */ + else + isteps = 16; + rsteps = 64; + errthr = 0.6; + if (doupdate) + mxits = 1; + else + mxits = 3; + break; + case 1: /* High */ + if (verify != 2 && doprofile && !doupdate) + isteps = 40; /* Use more steps if we're creating a profile */ + else + isteps = 20; + rsteps = 96; + errthr = 0.4; + if (doupdate) + mxits = 1; + else + mxits = 4; + break; + case 2: /* Ultra */ + if (verify != 2 && doprofile && !doupdate) + isteps = 48; /* Use more steps if we're creating a profile */ + else + isteps = 24; + rsteps = 128; + errthr = 0.25; + if (doupdate) + mxits = 1; + else + mxits = 5; + break; + } + + /* Set native white target flag in calx so that other things can play the game.. */ + if (wpx == 0.0 && wpy == 0.0 && temp == 0.0 && tbright == 0.0) + x.nat = 1; + else + x.nat = 0; + + /* Say something about what we're doing */ + if (verb) { + if (dtype > 0) + printf("Display type is '%c'\n",dtype); + + if (doupdate) { + if (x.nat) + printf("Target white = native white point & brightness\n"); + else + printf("Target white = XYZ %f %f %f\n", + x.twh[0], x.twh[1], x.twh[2]); + } else { + if (wpx > 0.0 || wpy > 0.0) + printf("Target white = xy %f %f\n",wpx,wpy); + else if (temp > 0.0) { + if (planckian) + printf("Target white = %f degrees kelvin Planckian (black body) spectrum\n",temp); + else + printf("Target white = %f degrees kelvin Daylight spectrum\n",temp); + } else + printf("Target white = native white point\n"); + + if (tbright > 0.0) + printf("Target white brightness = %f cd/m^2\n",tbright); + else + printf("Target white brightness = native brightness\n"); + if (bkbright > 0.0) + printf("Target black brightness = %f cd/m^2\n",bkbright); + else + printf("Target black brightness = native brightness\n"); + } + + switch(x.gammat) { + case gt_power: + if (egamma > 0.0) + printf("Target effective gamma = %f\n",egamma); + else + printf("Target advertised gamma = %f\n",gamma); + break; + case gt_Lab: + printf("Target gamma = L* curve\n"); + break; + case gt_sRGB: + printf("Target gamma = sRGB curve\n"); + break; + case gt_Rec709: + printf("Target gamma = REC 709 curve\n"); + break; + case gt_SMPTE240M: + printf("Target gamma = SMPTE 240M curve\n"); + break; + default: + error("Unknown gamma type"); + } + } + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* Go through the procedure of adjusting monitor controls */ + if (docontrols) { + int rgbch = 0; /* Got RBG Yxy ? */ + double rgbXYZ[3][3]; /* The RGB XYZ */ + + /* Make sure drift comp. is off for interactive adjustment */ + dr->change_drift_comp(dr, 0, 0); + + /* Until the user is done */ + printf("\nDisplay adjustment menu:"); + for (;;) { + int c; + + /* Print the menue of adjustments */ + printf("\nPress 1 .. 7\n"); + printf("1) Black level (CRT: Offset/Brightness)\n"); + printf("2) White point (Color temperature, R,G,B, Gain/Contrast)\n"); + printf("3) White level (CRT: Gain/Contrast, LCD: Brightness/Backlight)\n"); + printf("4) Black point (R,G,B, Offset/Brightness)\n"); + printf("5) Check all\n"); + printf("6) Measure and set ambient for viewing condition adjustment\n"); + printf("7) Continue on to calibration\n"); + printf("8) Exit\n"); + + empty_con_chars(); + c = next_con_char(); + + /* Black level adjustment */ + /* Due to the possibility of the channel offsets not being even, */ + /* we use the largest of the XYZ values after they have been */ + /* scaled to be even acording to the white XYZ balance. */ + /* It's safer to set the black level a bit low, and then the */ + /* calibration curves can bump the low ones up. */ + if (c == '1') { + col tcols[3] = { /* Base set of test colors */ + { 0.0, 0.0, 0.0 }, + { 0.5, 0.5, 0.5 }, /* And 1% values */ + { 1.0, 1.0, 1.0 } + }; + int ff; + double mgamma, tar1, dev1; + + printf("Doing some initial measurements\n"); + /* Do an initial set of readings to set 1% output mark */ + if ((rv = dr->read(dr, tcols, 3, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + + if (verb) { + printf("Black = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0], + tcols[0].XYZ[1], tcols[0].XYZ[2]); + printf("Grey = XYZ %6.2f %6.2f %6.2f\n",tcols[1].XYZ[0], + tcols[1].XYZ[1], tcols[1].XYZ[2]); + printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[2].XYZ[0], + tcols[2].XYZ[1], tcols[2].XYZ[2]); + } + + /* Advertised Gamma - Gross curve shape */ + mgamma = pop_gamma(tcols[0].XYZ[1], tcols[1].XYZ[1], tcols[2].XYZ[1]); + + dev1 = pow(0.01, 1.0/mgamma); +//printf("~1 device level for 1%% output = %f\n",dev1); + tcols[1].r = tcols[1].g = tcols[1].b = dev1; + tar1 = 0.01 * tcols[2].XYZ[1]; + + printf("\nAdjust CRT brightness to get target level. Press space when done.\n"); + printf(" Target %.2f\n",tar1); + for (ff = 0;; ff ^= 1) { + double dir; /* Direction to adjust brightness */ + double sv[3], val1; /* Scaled values */ + if ((rv = dr->read(dr, tcols+1, 1, 0, 0, 1, ' ',instClamp)) != 0) { + if (rv == 4) + break; /* User is done with this adjustment */ + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + /* Scale 1% values by ratio of Y to white XYZ */ + sv[0] = tcols[1].XYZ[0] * tcols[2].XYZ[1]/tcols[2].XYZ[0]; + sv[1] = tcols[1].XYZ[1]; + sv[2] = tcols[1].XYZ[2] * tcols[2].XYZ[1]/tcols[2].XYZ[2]; +//printf("~1 scaled readings = %f %f %f\n",sv[0],sv[1],sv[2]); + val1 = sv[1]; + if (sv[0] > val1) + val1 = sv[0]; + if (sv[2] > val1) + val1 = sv[2]; + dir = tar1 - val1; + if (fabs(dir) < 0.01) + dir = 0.0; + printf("%c%c Current %.2f %c", + cr_char, + ff == 0 ? '/' : '\\', + val1, + dir < 0.0 ? '-' : dir > 0.0 ? '+' : '='); + fflush(stdout); + } + printf("\n"); + + /* White point adjustment */ + } else if (c == '2') { + int nat = 0; /* NZ if using native white as target */ + col tcols[1] = { /* Base set of test colors */ + { 1.0, 1.0, 1.0 } + }; + int ff; + double tYxy[3]; /* Target white chromaticities */ + icmXYZNumber tXYZ; /* Target white as XYZ */ + double tLab[3]; /* Target white as Lab or UCS */ + double tarw; /* Target brightness */ + double Lab[3]; /* Last measured point Lab or UCS */ + double ct = 0.0, ct_de; /* Color temperature & delta E to white locus */ + + printf("Doing some initial measurements\n"); + + if (rgbch == 0) { /* Figure the RGB chromaticities */ + col ccols[3] = { + { 1.0, 0.0, 0.0 }, + { 0.0, 1.0, 0.0 }, + { 0.0, 0.0, 1.0 } + }; + if ((rv = dr->read(dr, ccols, 3, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("Red = XYZ %6.2f %6.2f %6.2f\n",ccols[0].XYZ[0], + ccols[0].XYZ[1], ccols[0].XYZ[2]); + printf("Green = XYZ %6.2f %6.2f %6.2f\n",ccols[1].XYZ[0], + ccols[1].XYZ[1], ccols[1].XYZ[2]); + printf("Blue = XYZ %6.2f %6.2f %6.2f\n",ccols[2].XYZ[0], + ccols[2].XYZ[1], ccols[2].XYZ[2]); + } + for (i = 0; i < 3; i++) + icmAry2Ary(rgbXYZ[i], ccols[i].XYZ); + rgbch = 1; + } + /* Do an initial set of readings to set full output mark */ + if ((rv = dr->read(dr, tcols, 1, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0], + tcols[0].XYZ[1], tcols[0].XYZ[2]); + } + + /* Figure out the target white chromaticity */ + if (wpx > 0.0 || wpy > 0.0) { /* xy coordinates */ + tYxy[0] = 1.0; + tYxy[1] = wpx; + tYxy[2] = wpy; + } else if (temp > 0.0) { /* Daylight color temperature */ + double XYZ[3]; + if (planckian) + rv = icx_ill_sp2XYZ(XYZ, icxOT_default, NULL, icxIT_Ptemp, temp, NULL); + else + rv = icx_ill_sp2XYZ(XYZ, icxOT_default, NULL, icxIT_Dtemp, temp, NULL); + if (rv != 0) + error("Failed to compute XYZ of target color temperature %f\n",temp); + icmXYZ2Yxy(tYxy, XYZ); + } else { /* Native white */ + icmXYZ2Yxy(tYxy, tcols[0].XYZ); + nat = 1; + } + + /* Figure out the target white brightness */ + /* Note we're not taking the device gamut into account here */ + if (tbright > 0.0) { /* Given brightness */ + tarw = tbright; +//printf("~1 set tarw %f from tbright\n",tarw); + } else { /* Native/maximum brightness */ + tarw = tcols[0].XYZ[1]; +//printf("~1 set tarw %f from tcols[1]\n",tarw); + } + + if (!nat) { /* Target is a specified white */ + printf("\nAdjust R,G & B gain to get target x,y. Press space when done.\n"); + printf(" Target Br %.2f, x %.4f , y %.4f \n", + tarw, tYxy[1],tYxy[2]); + + } else { /* Target is native white */ + printf("\nAdjust R,G & B gain to desired white point. Press space when done.\n"); + /* Compute the CT and delta E to white locus of target */ + ct = comp_ct(&ct_de, NULL, planckian, dovct, obType, tcols[0].XYZ); + printf(" Initial Br %.2f, x %.4f , y %.4f , %c%cT %4.0fK DE 2K %4.1f\n", + tarw, tYxy[1],tYxy[2], + dovct ? 'V' : 'C', planckian ? 'C' : 'D', ct,ct_de); + } + for (ff = 0;; ff ^= 1) { + double dir; /* Direction to adjust brightness */ + double Yxy[3]; /* Yxy of current reading */ + double rgbdir[3]; /* Direction to adjust RGB */ + double rgbxdir[3]; /* Biggest to move */ + double bdir, terr; + int bx = 0; + + if ((rv = dr->read(dr, tcols, 1, 0, 0, 1, ' ', instClamp)) != 0) { + if (rv == 4) + break; /* User is done with this adjustment */ + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + dir = tarw - tcols[0].XYZ[1]; + if (fabs(dir) < 0.01) + dir = 0.0; + + icmXYZ2Yxy(Yxy, tcols[0].XYZ); + + if (!nat) { /* Target is a specified white */ + /* Compute values we need for delta E and RGB direction */ + icmYxy2XYZ(tLab, tYxy); + tLab[0] /= tLab[1]; + tLab[2] /= tLab[1]; + tLab[1] /= tLab[1]; + icmAry2XYZ(tXYZ, tLab); /* Lab white reference */ + icmXYZ2Lab(&tXYZ, tLab, tLab); /* Target Lab */ + + icmAry2Ary(Lab, tcols[0].XYZ); + Lab[0] /= Lab[1]; + Lab[2] /= Lab[1]; + Lab[1] /= Lab[1]; + icmXYZ2Lab(&tXYZ, Lab, Lab); /* Current Lab */ + + } else { /* Target is native white */ + double lxyz[3]; /* Locus XYZ */ + ct = comp_ct(&ct_de, lxyz, planckian, dovct, obType, tcols[0].XYZ); + + icmXYZ2Yxy(tYxy, lxyz); + /* lxyz is already normalised */ + icmAry2XYZ(tXYZ, lxyz); /* Lab white reference */ + if (dovct) + icmXYZ2Lab(&tXYZ, tLab, lxyz); /* Target Lab */ + else + icmXYZ21960UCS(tLab, lxyz); /* Target UCS */ + + icmAry2Ary(Lab, tcols[0].XYZ); + Lab[0] /= Lab[1]; + Lab[2] /= Lab[1]; + Lab[1] /= Lab[1]; + if (dovct) + icmXYZ2Lab(&tXYZ, Lab, Lab); /* Current Lab */ + else + icmXYZ21960UCS(Lab, Lab); /* Current UCS */ + } + + /* Compute dot products */ + bdir = 0.0; + for (i = 0; i < 3; i++) { + double rgbLab[3]; + + if (dovct) + icmXYZ2Lab(&tXYZ, rgbLab, rgbXYZ[i]); + else + icmXYZ21960UCS(rgbLab, rgbXYZ[i]); + rgbdir[i] = (tLab[1] - Lab[1]) * (rgbLab[1] - Lab[1]) + + (tLab[2] - Lab[2]) * (rgbLab[2] - Lab[2]); + rgbxdir[i] = 0.0; + if (fabs(rgbdir[i]) > fabs(bdir)) { + bdir = rgbdir[i]; + bx = i; + } + } + + /* See how close to the target we are */ + terr = sqrt((tLab[1] - Lab[1]) * (tLab[1] - Lab[1]) + + (tLab[2] - Lab[2]) * (tLab[2] - Lab[2])); + if (terr < 0.1) + rgbdir[0] = rgbdir[1] = rgbdir[2] = 0.0; + rgbxdir[bx] = rgbdir[bx]; + + + if (!nat) { + printf("%c%c Current Br %.2f, x %.4f%c, y %.4f%c DE %4.1f R%c%c G%c%c B%c%c ", + cr_char, + ff == 0 ? '/' : '\\', + tcols[0].XYZ[1], + Yxy[1], + Yxy[1] > tYxy[1] ? '-' : Yxy[1] < tYxy[1] ? '+' : '=', + Yxy[2], + Yxy[2] > tYxy[2] ? '-' : Yxy[2] < tYxy[2] ? '+' : '=', + icmCIE2K(tLab, Lab), + rgbdir[0] < 0.0 ? '-' : rgbdir[0] > 0.0 ? '+' : '=', + rgbxdir[0] < 0.0 ? '-' : rgbxdir[0] > 0.0 ? '+' : ' ', + rgbdir[1] < 0.0 ? '-' : rgbdir[1] > 0.0 ? '+' : '=', + rgbxdir[1] < 0.0 ? '-' : rgbxdir[1] > 0.0 ? '+' : ' ', + rgbdir[2] < 0.0 ? '-' : rgbdir[2] > 0.0 ? '+' : '=', + rgbxdir[2] < 0.0 ? '-' : rgbxdir[2] > 0.0 ? '+' : ' '); + } else { + printf("%c%c Current Br %.2f, x %.4f%c, y %.4f%c %c%cT %4.0fK DE 2K %4.1f R%c%c G%c%c B%c%c ", + cr_char, + ff == 0 ? '/' : '\\', + tcols[0].XYZ[1], + Yxy[1], + Yxy[1] > tYxy[1] ? '-': Yxy[1] < tYxy[1] ? '+' : '=', + Yxy[2], + Yxy[2] > tYxy[2] ? '-': Yxy[2] < tYxy[2] ? '+' : '=', + dovct ? 'V' : 'C', planckian ? 'C' : 'D', ct,ct_de, + rgbdir[0] < 0.0 ? '-' : rgbdir[0] > 0.0 ? '+' : '=', + rgbxdir[0] < 0.0 ? '-' : rgbxdir[0] > 0.0 ? '+' : ' ', + rgbdir[1] < 0.0 ? '-' : rgbdir[1] > 0.0 ? '+' : '=', + rgbxdir[1] < 0.0 ? '-' : rgbxdir[1] > 0.0 ? '+' : ' ', + rgbdir[2] < 0.0 ? '-' : rgbdir[2] > 0.0 ? '+' : '=', + rgbxdir[2] < 0.0 ? '-' : rgbxdir[2] > 0.0 ? '+' : ' '); + } + fflush(stdout); + } + printf("\n"); + + /* White level adjustment */ + } else if (c == '3') { + col tcols[1] = { /* Base set of test colors */ + { 1.0, 1.0, 1.0 } + }; + int ff; + double tarw; + + printf("Doing some initial measurements\n"); + /* Do an initial set of readings to set full output mark */ + if ((rv = dr->read(dr, tcols, 1, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0], + tcols[0].XYZ[1], tcols[0].XYZ[2]); + } + + /* Figure out the target white brightness */ + /* Note we're not taking the device gamut into account here */ + if (tbright > 0.0) /* Given brightness */ + tarw = tbright; + else /* Native/maximum brightness */ + tarw = tcols[0].XYZ[1]; + + if (tbright > 0.0) { + printf("\nAdjust CRT Contrast or LCD Brightness to get target level. Press space when done.\n"); + printf(" Target %.2f\n", tarw); + } else { + printf("\nAdjust CRT Contrast or LCD Brightness to desired level. Press space when done.\n"); + printf(" Initial %.2f\n", tarw); + } + for (ff = 0;; ff ^= 1) { + double dir; /* Direction to adjust brightness */ + + if ((rv = dr->read(dr, tcols, 1, 0, 0, 1, ' ', instClamp)) != 0) { + if (rv == 4) + break; /* User is done with this adjustment */ + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + dir = tarw - tcols[0].XYZ[1]; + if (fabs(dir) < 0.01) + dir = 0.0; + + if (tbright > 0.0) + printf("%c%c Current %.2f %c", + cr_char, + ff == 0 ? '/' : '\\', + tcols[0].XYZ[1], + dir < 0.0 ? '-' : dir > 0.0 ? '+' : '='); + else + printf("%c%c Current %.2f ", + cr_char, + ff == 0 ? '/' : '\\', + tcols[0].XYZ[1]); + fflush(stdout); + } + printf("\n"); + + /* Black point adjustment */ + } else if (c == '4') { + col tcols[3] = { /* Base set of test colors */ + { 0.0, 0.0, 0.0 }, + { 0.5, 0.5, 0.5 }, /* And 1% values */ + { 1.0, 1.0, 1.0 } + }; + int ff; + double tYxy[3]; /* Target white chromaticities */ + icmXYZNumber tXYZ; /* Target white as XYZ */ + double tLab[3]; /* Target white as Lab or UCS */ + double mgamma, tar1, dev1; + double Lab[3]; /* Last measured point Lab */ + + printf("Doing some initial measurements\n"); + + if (rgbch == 0) { /* Figure the RGB chromaticities */ + col ccols[3] = { + { 1.0, 0.0, 0.0 }, + { 0.0, 1.0, 0.0 }, + { 0.0, 0.0, 1.0 } + }; + if ((rv = dr->read(dr, ccols, 3, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("Red = XYZ %6.2f %6.2f %6.2f\n",ccols[0].XYZ[0], + ccols[0].XYZ[1], ccols[0].XYZ[2]); + printf("Green = XYZ %6.2f %6.2f %6.2f\n",ccols[1].XYZ[0], + ccols[1].XYZ[1], ccols[1].XYZ[2]); + printf("Blue = XYZ %6.2f %6.2f %6.2f\n",ccols[2].XYZ[0], + ccols[2].XYZ[1], ccols[2].XYZ[2]); + } + for (i = 0; i < 3; i++) + icmAry2Ary(rgbXYZ[i], ccols[i].XYZ); + rgbch = 1; + } + /* Do an initial set of readings to set 1% output mark */ + if ((rv = dr->read(dr, tcols, 3, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("Black = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0], + tcols[0].XYZ[1], tcols[0].XYZ[2]); + printf("Grey = XYZ %6.2f %6.2f %6.2f\n",tcols[1].XYZ[0], + tcols[1].XYZ[1], tcols[1].XYZ[2]); + printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[2].XYZ[0], + tcols[2].XYZ[1], tcols[2].XYZ[2]); + } + + /* Advertised Gamma - Gross curve shape */ + mgamma = pop_gamma(tcols[0].XYZ[1], tcols[1].XYZ[1], tcols[2].XYZ[1]); + + dev1 = pow(0.01, 1.0/mgamma); + tcols[1].r = tcols[1].g = tcols[1].b = dev1; + tar1 = 0.01 * tcols[2].XYZ[1]; + + /* Figure out the target white chromaticity */ + if (wpx > 0.0 || wpy > 0.0) { /* xy coordinates */ + tYxy[0] = 1.0; + tYxy[1] = wpx; + tYxy[2] = wpy; + + } else if (temp > 0.0) { /* Daylight color temperature */ + double XYZ[3]; + if (planckian) + rv = icx_ill_sp2XYZ(XYZ, icxOT_default, NULL, icxIT_Ptemp, temp, NULL); + else + rv = icx_ill_sp2XYZ(XYZ, icxOT_default, NULL, icxIT_Dtemp, temp, NULL); + if (rv != 0) + error("Failed to compute XYZ of target color temperature %f\n",temp); + icmXYZ2Yxy(tYxy, XYZ); + } else { /* Native white */ + icmXYZ2Yxy(tYxy, tcols[2].XYZ); + } + + printf("\nAdjust R,G & B offsets to get target x,y. Press space when done.\n"); + printf(" Target Br %.2f, x %.4f , y %.4f \n", tar1, tYxy[1],tYxy[2]); + for (ff = 0;; ff ^= 1) { + double dir; /* Direction to adjust brightness */ + double sv[3], val1; /* Scaled values */ + double Yxy[3]; /* Yxy of current reading */ + double rgbdir[3]; /* Direction to adjust RGB */ + double rgbxdir[3]; /* Biggest to move */ + double bdir, terr; + int bx = 0; + + if ((rv = dr->read(dr, tcols+1, 1, 0, 0, 1, ' ', instClamp)) != 0) { + if (rv == 4) + break; /* User is done with this adjustment */ + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + + /* Scale 1% values by ratio of Y to white XYZ */ + sv[0] = tcols[1].XYZ[0] * tcols[2].XYZ[1]/tcols[2].XYZ[0]; + sv[1] = tcols[1].XYZ[1]; + sv[2] = tcols[1].XYZ[2] * tcols[2].XYZ[1]/tcols[2].XYZ[2]; + val1 = sv[1]; + if (sv[0] > val1) + val1 = sv[0]; + if (sv[2] > val1) + val1 = sv[2]; + + /* Compute 1% direction */ + dir = tar1 - val1; + if (fabs(dir) < 0.01) + dir = 0.0; + + /* Compute numbers for black point error and direction */ + icmYxy2XYZ(tLab, tYxy); + tLab[0] /= tLab[1]; + tLab[2] /= tLab[1]; + tLab[1] /= tLab[1]; + icmAry2XYZ(tXYZ, tLab); /* Lab white reference */ + icmXYZ2Lab(&tXYZ, tLab, tLab); + + icmXYZ2Yxy(Yxy, tcols[1].XYZ); + icmAry2Ary(Lab, tcols[1].XYZ); + Lab[0] /= Lab[1]; + Lab[2] /= Lab[1]; + Lab[1] /= Lab[1]; + icmXYZ2Lab(&tXYZ, Lab, Lab); + + /* Compute dot products */ + bdir = 0.0; + for (i = 0; i < 3; i++) { + double rgbLab[3]; + + icmXYZ2Lab(&tXYZ, rgbLab, rgbXYZ[i]); + rgbdir[i] = (tLab[1] - Lab[1]) * (rgbLab[1] - Lab[1]) + + (tLab[2] - Lab[2]) * (rgbLab[2] - Lab[2]); + rgbxdir[i] = 0.0; + if (fabs(rgbdir[i]) > fabs(bdir)) { + bdir = rgbdir[i]; + bx = i; + } + } + + /* See how close to the target we are */ + terr = sqrt((tLab[1] - Lab[1]) * (tLab[1] - Lab[1]) + + (tLab[2] - Lab[2]) * (tLab[2] - Lab[2])); + if (terr < 0.1) + rgbdir[0] = rgbdir[1] = rgbdir[2] = 0.0; + rgbxdir[bx] = rgbdir[bx]; + + printf("%c%c Current Br %.2f, x %.4f%c, y %.4f%c DE %4.1f R%c%c G%c%c B%c%c ", + cr_char, + ff == 0 ? '/' : '\\', + val1, + Yxy[1], + Yxy[1] > tYxy[1] ? '-': Yxy[1] < tYxy[1] ? '+' : '=', + Yxy[2], + Yxy[2] > tYxy[2] ? '-': Yxy[2] < tYxy[2] ? '+' : '=', + icmCIE2K(tLab, Lab), + rgbdir[0] < 0.0 ? '-' : rgbdir[0] > 0.0 ? '+' : '=', + rgbxdir[0] < 0.0 ? '-' : rgbxdir[0] > 0.0 ? '+' : ' ', + rgbdir[1] < 0.0 ? '-' : rgbdir[1] > 0.0 ? '+' : '=', + rgbxdir[1] < 0.0 ? '-' : rgbxdir[1] > 0.0 ? '+' : ' ', + rgbdir[2] < 0.0 ? '-' : rgbdir[2] > 0.0 ? '+' : '=', + rgbxdir[2] < 0.0 ? '-' : rgbxdir[2] > 0.0 ? '+' : ' '); + fflush(stdout); + } + printf("\n"); + + /* Report on how well we current meet the targets */ + } else if (c == '5') { + int nat = 0; /* NZ if using native white as target */ + col tcols[4] = { /* Set of test colors */ + { 0.0, 0.0, 0.0 }, + { 0.5, 0.5, 0.5 }, + { 1.0, 1.0, 1.0 }, + { 0.0, 0.0, 0.0 } /* 1% test value */ + }; + double tYxy[3]; /* Target white chromaticities */ + double sv[3], val1; /* Scaled values */ + double mgamma, tarw, tar1, dev1, tarh; + double gooff; /* Aproximate output offset needed */ + icmXYZNumber tXYZ; + double tLab[3], wYxy[3], wLab[3], bYxy[3], bLab[3]; + double ct, ct_de; /* Color temperature & delta E to white locus */ + + printf("Doing check measurements\n"); + + /* Do an initial set of readings to set 1% output mark */ + if ((rv = dr->read(dr, tcols, 3, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("Black = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0], + tcols[0].XYZ[1], tcols[0].XYZ[2]); + printf("Grey = XYZ %6.2f %6.2f %6.2f\n",tcols[1].XYZ[0], + tcols[1].XYZ[1], tcols[1].XYZ[2]); + printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[2].XYZ[0], + tcols[2].XYZ[1], tcols[2].XYZ[2]); + } + + /* Approximate Gamma - use the gross curve shape for robustness */ + mgamma = pop_gamma(tcols[0].XYZ[1], tcols[1].XYZ[1], tcols[2].XYZ[1]); + + dev1 = pow(0.01, 1.0/mgamma); + tcols[3].r = tcols[3].g = tcols[3].b = dev1; + tar1 = 0.01 * tcols[2].XYZ[1]; + + /* Read the 1% value */ + if ((rv = dr->read(dr, tcols+3, 1, 0, 0, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + if (verb) { + printf("1%% = XYZ %6.2f %6.2f %6.2f\n",tcols[3].XYZ[0], + tcols[3].XYZ[1], tcols[3].XYZ[2]); + } + + /* Scale 1% values by ratio of Y to white XYZ */ + /* (Note we're assuming -k1 here, which may not be true...) */ + sv[0] = tcols[3].XYZ[0] * tcols[2].XYZ[1]/tcols[2].XYZ[0]; + sv[1] = tcols[3].XYZ[1]; + sv[2] = tcols[3].XYZ[2] * tcols[2].XYZ[1]/tcols[2].XYZ[2]; + val1 = sv[1]; + if (sv[0] > val1) + val1 = sv[0]; + if (sv[2] > val1) + val1 = sv[2]; + + /* Figure out the target white brightness */ + /* Note we're not taking the device gamut into account here */ + if (tbright > 0.0) /* Given brightness */ + tarw = tbright; + else /* Native/maximum brightness */ + tarw = tcols[2].XYZ[1]; + + /* Figure out the target white chromaticity */ + if (wpx > 0.0 || wpy > 0.0) { /* xy coordinates */ + tYxy[0] = 1.0; + tYxy[1] = wpx; + tYxy[2] = wpy; + + } else if (temp > 0.0) { /* Daylight color temperature */ + double XYZ[3]; + if (planckian) + rv = icx_ill_sp2XYZ(XYZ, icxOT_default, NULL, icxIT_Ptemp, temp, NULL); + else + rv = icx_ill_sp2XYZ(XYZ, icxOT_default, NULL, icxIT_Dtemp, temp, NULL); + if (rv != 0) + error("Failed to compute XYZ of target color temperature %f\n",temp); + icmXYZ2Yxy(tYxy, XYZ); + } else { /* Native white */ + icmXYZ2Yxy(tYxy, tcols[2].XYZ); + nat = 1; + } + + /* Figure out the target 50% device output value */ + gooff = tcols[0].XYZ[1]/tcols[2].XYZ[1]; /* Aprox. normed black output offset */ + + /* Use tech_gamma() to do the hard work */ + tarh = tech_gamma(&x, NULL, NULL, NULL, egamma, gamma, gooff); + + /* Convert from Y fraction to absolute Y */ + tarh = tarh * tcols[2].XYZ[1]; + + /* Compute various white point values */ + icmYxy2XYZ(tLab, tYxy); + tLab[0] /= tLab[1]; + tLab[2] /= tLab[1]; + tLab[1] /= tLab[1]; + icmAry2XYZ(tXYZ, tLab); + icmXYZ2Lab(&tXYZ, tLab, tLab); + + icmXYZ2Yxy(wYxy, tcols[2].XYZ); + icmAry2Ary(wLab, tcols[2].XYZ); + wLab[0] /= wLab[1]; + wLab[2] /= wLab[1]; + wLab[1] /= wLab[1]; + icmXYZ2Lab(&tXYZ, wLab, wLab); + + icmXYZ2Yxy(bYxy, tcols[3].XYZ); + icmAry2Ary(bLab, tcols[3].XYZ); + bLab[0] /= bLab[1]; + bLab[2] /= bLab[1]; + bLab[1] /= bLab[1]; + icmXYZ2Lab(&tXYZ, bLab, bLab); + + /* And color temperature */ + ct = comp_ct(&ct_de, NULL, planckian, dovct, obType, tcols[2].XYZ); + + printf("\n"); + + if (tbright > 0.0) /* Given brightness */ + printf(" Target Brightness = %.2f, Current = %5.2f, error = % .1f%%\n", + tarw, tcols[2].XYZ[1], + 100.0 * (tcols[2].XYZ[1] - tarw)/tarw); + else + printf(" Current Brightness = %.2f\n", tcols[2].XYZ[1]); + + printf(" Target 50%% Level = %.2f, Current = %5.2f, error = % .1f%%\n", + tarh, tcols[1].XYZ[1], + 100.0 * (tcols[1].XYZ[1] - tarh)/tarw); + + printf(" Target Near Black = %5.2f, Current = %5.2f, error = % .1f%%\n", + tar1, val1, + 100.0 * (val1 - tar1)/tarw); + + if (!nat) + printf(" Target white = x %.4f, y %.4f, Current = x %.4f, y %.4f, error = %5.2f DE\n", + tYxy[1], tYxy[2], wYxy[1], wYxy[2], icmCIE2K(tLab, wLab)); + else + printf(" Current white = x %.4f, y %.4f, %c%cT %4.0fK DE 2K %4.1f\n", + wYxy[1], wYxy[2], dovct ? 'V' : 'C', planckian ? 'C' : 'D', ct,ct_de); + + printf(" Target black = x %.4f, y %.4f, Current = x %.4f, y %.4f, error = %5.2f DE\n", + tYxy[1], tYxy[2], bYxy[1], bYxy[2], icmCIE2K(tLab, bLab)); + + + /* Measure and set ambient for viewing condition adjustment */ + } else if (c == '6') { + if ((rv = dr->ambient(dr, &ambient, 1)) != 0) { + if (rv == 8) { + printf("Instrument doesn't have an ambient reading capability\n"); + } else { + dr->del(dr); + error("ambient measure failed with '%s'\n",disprd_err(rv)); + } + } else { + printf("Measured ambient level = %.1f Lux\n",ambient * 3.141592654); + } + + } else if (c == '7') { + if (!verb) { /* Tell user command has been accepted */ + if (verify == 2) + printf("Commencing device verification\n"); + else + printf("Commencing device calibration\n"); + } + break; + } else if (c == '8' || c == 0x03 || c == 0x1b) { + printf("Exiting\n"); + dr->del(dr); + exit(0); + } + } + + /* Make sure drift comp. is set to the command line options */ + dr->change_drift_comp(dr, bdrift, wdrift); + } + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + /* Take a small number of readings, and compute basic */ + /* informations such as black & white, white target, */ + /* aproximate matrix based display forward and reverse model. */ + /* If bkcorrect is auto, determine a level. */ + + /* Read the base test set */ + { + icmXYZNumber mrd; /* Number for matrix */ + icmXYZNumber mgn; + icmXYZNumber mbl; + icmXYZNumber mwh; + ramdac *or = NULL; + + col base[6] = { /* Base set of test colors */ + { 0.0, 0.0, 0.0 }, /* 0 - Black */ + { 1.0, 0.0, 0.0 }, /* 1 - Red */ + { 0.0, 1.0, 0.0 }, /* 2 - Green */ + { 0.0, 0.0, 1.0 }, /* 3 - Blue */ + { 1.0, 1.0, 1.0 }, /* 4 - White */ + { 0.0, 0.0, 0.0 } /* 5 - Black */ + }; + + if (verb) { + if (verify == 2) + printf("Commencing device verification\n"); + else + printf("Commencing device calibration\n"); + } + + /* Switch to native for this, so the black calc is realistic. */ + /* (Should we really get black aim from previous .cal though ???) */ + if (verify == 2) { + if ((or = dr->dw->get_ramdac(dr->dw)) != NULL) { + ramdac *r; + if (verb) printf("Switching to native response for base measurements\n"); + r = or->clone(or); + r->setlin(r); + dr->dw->set_ramdac(dr->dw, r, 0); + r->del(r); + } + } + + /* Read the patches without clamping */ + if ((rv = dr->read(dr, base, 6, 1, 6, 1, 0, instNoClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + + /* Restore the cal we're verifying */ + if (verify == 2 && or != NULL) { + if (verb) printf("Switching back to calibration being verified\n"); + dr->dw->set_ramdac(dr->dw, or, 0); + or->del(or); + } + + if (base[0].XYZ_v == 0) { + dr->del(dr); + error("Failed to get an XYZ value from the instrument!\n"); + } + + /* Average black relative from 2 readings */ + x.bk[0] = 0.5 * (base[0].XYZ[0] + base[5].XYZ[0]); + x.bk[1] = 0.5 * (base[0].XYZ[1] + base[5].XYZ[1]); + x.bk[2] = 0.5 * (base[0].XYZ[2] + base[5].XYZ[2]); + icmClamp3(x.bk, x.bk); /* And clamp them */ + for (i = 1; i < 5; i++) + icmClamp3(base[i].XYZ, base[i].XYZ); + + /* Copy other readings into place */ + dispLum = base[4].XYZ[1]; /* White Y */ + icmAry2Ary(x.wh, base[4].XYZ); + icmAry2XYZ(x.twN, x.wh); /* Use this as Lab reference white until we establish target */ + + icmAry2XYZ(mrd, base[1].XYZ); + icmAry2XYZ(mgn, base[2].XYZ); + icmAry2XYZ(mbl, base[3].XYZ); + icmAry2XYZ(mwh, base[4].XYZ); + + if (verb) { + printf("Black = XYZ %6.2f %6.2f %6.2f\n",x.bk[0],x.bk[1],x.bk[2]); + printf("Red = XYZ %6.2f %6.2f %6.2f\n",base[1].XYZ[0], base[1].XYZ[1], base[1].XYZ[2]); + printf("Green = XYZ %6.2f %6.2f %6.2f\n",base[2].XYZ[0], base[2].XYZ[1], base[2].XYZ[2]); + printf("Blue = XYZ %6.2f %6.2f %6.2f\n",base[3].XYZ[0], base[3].XYZ[1], base[3].XYZ[2]); + printf("White = XYZ %6.2f %6.2f %6.2f\n",base[4].XYZ[0], base[4].XYZ[1], base[4].XYZ[2]); + } + + /* Setup forward matrix */ + if (icmRGBprim2matrix(mwh, mrd, mgn, mbl, x.fm)) { + dr->del(dr); + error("Aprox. fwd matrix unexpectedly singular\n"); + } + +#ifdef DEBUG + if (verb) { + printf("Forward matrix is:\n"); + printf("%f %f %f\n", x.fm[0][0], x.fm[0][1], x.fm[0][2]); + printf("%f %f %f\n", x.fm[1][0], x.fm[1][1], x.fm[1][2]); + printf("%f %f %f\n", x.fm[2][0], x.fm[2][1], x.fm[2][2]); + } +#endif + + /* Compute bwd matrix */ + if (icmInverse3x3(x.bm, x.fm)) { + dr->del(dr); + error("Inverting aprox. fwd matrix failed"); + } + + /* Decide on the level of black correction. */ + if (bkcorrect < 0.0) { + double rat; + + /* rat is 0 for displays with a good black, */ + /* and 1 for displays with a bad black level. */ + /* (Not sure if this should be scaled by the white, */ + /* making it contrast ratio sensitive?) */ + rat = (x.bk[1] - 0.02)/(0.3 - 0.02); + if (rat < 0.0) + rat = 0.0; + else if (rat > 1.0) + rat = 1.0; + /* Make transition more perceptual */ + rat = sqrt(rat); + bkcorrect = 1.0 - rat; + if (verb) + printf("Automatic black point hue correction level = %1.2f\n", bkcorrect); + } + } + + /* Now do some more readings, to compute the basic per channel */ + /* transfer characteristics, and then a device model. */ + if (verify != 2 && !doupdate) { + col *cols; /* Read 4 x isteps patches from display */ + sxyz *asrgb[4]; /* samples for r, g, b & w */ + + if ((cols = (col *)malloc(isteps * 4 * sizeof(col))) == NULL) { + dr->del(dr); + error ("Malloc of array of readings failed"); + } + for (j = 0; j < 4; j++) { + if ((asrgb[j] = (sxyz *)malloc(isteps * sizeof(sxyz))) == NULL) { + free(cols); + dr->del(dr); + error ("Malloc of array of readings failed"); + } + } + + /* Set the device colors to read */ + for (i = 0; i < isteps; i++) { + double vv; + +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(i); +#endif + vv = i/(isteps - 1.0); + vv = pow(vv, MOD_DIST_POW); + for (j = 0; j < 4; j++) { + cols[i * 4 + j].r = cols[i * 4 + j].g = cols[i * 4 + j].b = 0.0; + if (j == 0) + cols[i * 4 + j].r = vv; + else if (j == 1) + cols[i * 4 + j].g = vv; + else if (j == 2) + cols[i * 4 + j].b = vv; + else + cols[i * 4 + j].r = cols[i * 4 + j].g = cols[i * 4 + j].b = vv; + } + } + + /* Read the patches */ + if ((rv = dr->read(dr, cols, isteps * 4, 1, isteps * 4, 1, 0, instClamp)) != 0) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + + /* Transfer readings to asrgb[] */ + for (i = 0; i < isteps; i++) { + double vv = cols[i * 4 + 0].r; + for (j = 0; j < 4; j++) { +//printf("~1 R = %f, G = %f, B = %f, XYZ = %f %f %f\n", +//cols[i * 4 + j].r, cols[i * 4 + j].g, cols[i * 4 + j].b, cols[i * 4 + j].XYZ[0], cols[i * 4 + j].XYZ[1], cols[i * 4 + j].XYZ[2]); + asrgb[j][i].v = vv; + asrgb[j][i].xyz[0] = cols[i * 4 + j].XYZ[0]; + asrgb[j][i].xyz[1] = cols[i * 4 + j].XYZ[1]; + asrgb[j][i].xyz[2] = cols[i * 4 + j].XYZ[2]; + } + } + + /* Convert RGB channel samples to curves */ + { + mcvco *sdv; /* Points used to create cvs[], RGB */ + double blrgb[3]; + double *op; /* Parameters to optimise */ + double *sa; /* Search area */ + double re; /* Residual error */ + + /* Transform measured black back to linearised RGB values */ + icmMulBy3x3(blrgb, x.bm, x.bk); +//printf("~1 model black should be %f %f %f\n", x.bk[0], x.bk[1], x.bk[2]); +//printf("~1 linearised RGB should be %f %f %f\n", blrgb[0], blrgb[1], blrgb[2]); + + if ((sdv = malloc(sizeof(mcvco) * isteps)) == NULL) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error ("Malloc of scattered data points failed"); + } + for (k = 0; k < 3; k++) { /* Create the model curves */ + for (i = 0; i < isteps; i++) { + sdv[i].p = asrgb[k][i].v; + sdv[i].v = ICMNORM3(asrgb[k][i].xyz); + sdv[i].w = 1.0; +//printf("~1 chan %d, entry %d, p = %f, v = %f from XYZ %f %f %f\n", +//k,i,x.sdv[k][i].p,x.sdv[k][i].v, asrgb[k][i].xyz[0], asrgb[k][i].xyz[1], asrgb[k][i].xyz[2]); + } + if ((x.dcvs[k] = new_mcv()) == NULL) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error("new_mcv x.dcvs[%d] failed",k); + } + x.dcvs[k]->fit(x.dcvs[k], 0, fitord, sdv, isteps, 5.0); + + /* Scale the whole curve so the output is scaled to 1.0 */ + x.dcvs[k]->force_scale(x.dcvs[k], 1.0); + + /* Force curves to produce this lrgb for 0.0 */ + x.dcvs[k]->force_0(x.dcvs[k], blrgb[k]); + } + free(sdv); + +#ifdef OPTIMIZE_MODEL + /* Setup list of reference points ready for optimisation */ + x.nrp = 4 * isteps; + if ((x.rp = (optref *)malloc(sizeof(optref) * x.nrp)) == NULL) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error ("Malloc of measurement reference points failed"); + } + for (k = 0; k < 4; k++) { + for (i = 0; i < isteps; i++) { + int ii = k * isteps + i; + double v[3]; + + v[0] = v[1] = v[2] = 0.0; + if (k == 0) + v[k] = asrgb[k][i].v; + else if (k == 1) + v[k] = asrgb[k][i].v; + else if (k == 2) + v[k] = asrgb[k][i].v; + else + v[0] = v[1] = v[2] = asrgb[k][i].v; + icmAry2Ary(x.rp[ii].dev, v); + icmXYZ2Lab(&x.twN, x.rp[ii].lab, asrgb[k][i].xyz); + if (k == 3) /* White */ + x.rp[ii].w = 0.5; + else + x.rp[ii].w = 0.16667; + } + } + + /* Get parameters and setup for optimisation */ + op = dev_get_params(&x); + if ((sa = malloc(x.np * sizeof(double))) == NULL) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error ("Malloc of scattered data points failed"); + } + + for (i = 0; i < x.np; i++) + sa[i] = 0.1; + + /* Do optimisation */ +#ifdef NEVER + if (powell(&re, x.np, op, sa, 1e-5, 3000, dev_opt_func, (void *)&x, NULL, NULL) != 0) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error ("Model powell failed, re = %f",re); + } +#else + if (conjgrad(&re, x.np, op, sa, 1e-5, 3000, + dev_opt_func, dev_dopt_func, (void *)&x, NULL, NULL) != 0) { + if (re > 1e-2) { + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + dr->del(dr); + error("Model conjgrad failed, residual error = %f",re); + } else + warning("Model conjgrad failed, residual error = %f",re); + } +#endif + + /* Put optimised parameters in place */ + dev_put_params(&x, op); + + free(x.rp); + x.rp = NULL; + x.nrp = 0; + free(x.dtin_iv); /* Free temporary arrays */ + x.dtin_iv = NULL; + free(sa); + free(op); +#endif /* OPTIMIZE_MODEL */ + } + +#ifdef DEBUG_PLOT + /* Plot the current calc curves */ + { + #define XRES 256 + double xx[XRES]; + double yy[3][XRES]; + double xyz[3]; + for (i = 0; i < XRES; i++) { + xx[i] = i/(XRES-1.0); + for (j = 0; j < 3; j++) + yy[j][i] = x.dcvs[j]->interp(x.dcvs[j], xx[i]); + } + printf("Channel curves\n"); + do_plot(xx,yy[0],yy[1],yy[2],XRES); + #undef XRES + } +#endif + + /* We're done with cols[] and asrgb[] */ + free(cols); free(asrgb[0]); free(asrgb[1]); free(asrgb[2]); free(asrgb[3]); + } + +#ifdef CHECK_MODEL + /* Check how well our fwd model agrees with the device */ + if (verify != 2) { + col set[3]; /* Variable to read up to 3 values from the display */ + int nn = 27; + double alab[3], mxyz[3], mlab[3]; /* Actual and model Lab */ + double mnerr; /* Maximum neutral error */ + double mnv; /* Value where maximum error is */ + double anerr; /* Average neutral error */ + + mnerr = anerr = 0.0; + /* !!! Should change this to single batch to work better with drift comp. !!! */ + for (i = 0; i < (nn + 3); i++) { + double vv, v[3]; + double de; + + if (i < nn) { + vv = i/(nn - 1.0); + vv = pow(vv, CHECK_DIST_POW); + v[0] = v[1] = v[2] = vv; + set[0].r = set[0].g = set[0].b = vv; + + } else { /* Do R, G, B */ + v[0] = v[1] = v[2] = 0.0; + v[i - nn] = 1.0; + set[0].r = v[0]; + set[0].g = v[1]; + set[0].b = v[2]; + } + + if ((rv = dr->read(dr, set, 1, i+1, nn+3, 1, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + icmXYZ2Lab(&x.twN, alab, set[0].XYZ); + + fwddev(&x, mxyz, v); + icmXYZ2Lab(&x.twN, mlab, mxyz); + + de = icmCIE2K(mlab, alab); + if (de > mnerr) { + mnerr = de; + mnv = vv; + } + anerr += de; + + printf("RGB %.3f %.3f %.3f -> XYZ %.2f %.2f %.2f, model %.2f %.2f %.2f\n", + set[0].r, set[0].g, set[0].b, + set[0].XYZ[0], set[0].XYZ[1], + set[0].XYZ[2], mxyz[0], mxyz[1], mxyz[2]); + + printf("RGB %.3f %.3f %.3f -> Lab %.2f %.2f %.2f, model %.2f %.2f %.2f, DE %f\n", + set[0].r, set[0].g, set[0].b, alab[0], alab[1], alab[2], mlab[0], mlab[1], mlab[2],de); + } + anerr /= (double)(nn+3); + printf("Model maximum error (@ %f) = %f deltaE\n",mnv, mnerr); + printf("Model average error = %f deltaE\n",anerr); + } +#endif /* CHECK_MODEL */ + + /* Figure out our calibration curve parameter targets */ + if (!doupdate) { + + /* Figure out the target white point */ + if (wpx > 0.0 || wpy > 0.0) { /* xy coordinates */ + double Yxy[3]; + Yxy[0] = 1.0; + Yxy[1] = wpx; + Yxy[2] = wpy; + icmYxy2XYZ(x.twh, Yxy); + + } else if (temp > 0.0) { /* Daylight color temperature */ + if (planckian) + rv = icx_ill_sp2XYZ(x.twh, icxOT_default, NULL, icxIT_Ptemp, temp, NULL); + else + rv = icx_ill_sp2XYZ(x.twh, icxOT_default, NULL, icxIT_Dtemp, temp, NULL); + if (rv != 0) + error("Failed to compute XYZ of target color temperature %f\n",temp); +//printf("~1 Raw target from temp %f XYZ = %f %f %f\n",temp,x.twh[0],x.twh[1],x.twh[2]); + } else { /* Native white */ + x.twh[0] = x.wh[0]/x.wh[1]; + x.twh[1] = x.wh[1]/x.wh[1]; + x.twh[2] = x.wh[2]/x.wh[1]; + } + x.nwh[0] = x.twh[0]; + x.nwh[1] = x.twh[1]; + x.nwh[2] = x.twh[2]; + + /* Convert it to absolute white target */ + if (tbright > 0.0) { /* Given brightness */ + x.twh[0] *= tbright; + x.twh[1] *= tbright; + x.twh[2] *= tbright; + } else { /* Native/maximum brightness */ + x.twh[0] *= x.wh[1]; + x.twh[1] *= x.wh[1]; + x.twh[2] *= x.wh[1]; + if (verb) + printf("Initial native brightness target = %f cd/m^2\n", x.twh[1]); + } + + /* Now make sure the target white will fit in gamut. */ + if (verify != 2 && + ((tbright > 0.0 && invlindev(&x, NULL, x.twh) > 0.0) /* Defined brightness and clips */ + || (tbright <= 0.0 && x.nat == 0))) { /* Max non-native white */ + double rgb[3]; + double scale = 0.5; + double sa = 0.1; + + if (powell(NULL, 1, &scale, &sa, 1e-7, 500, wp_opt_func, (void *)&x, NULL, NULL) != 0) + error ("WP scale powell failed"); + + x.twh[0] *= scale; + x.twh[1] *= scale; + x.twh[2] *= scale; + invdev(&x, rgb, x.twh); + if (verb) { + printf("Had to scale brightness from %f to %f to fit within gamut,\n",x.twh[1]/scale, x.twh[1]); + printf("corresponding to aprox. RGB %f %f %f\n",rgb[0],rgb[1],rgb[2]); + } + } + + if (verb) + printf("Target white value is XYZ %f %f %f\n",x.twh[0],x.twh[1],x.twh[2]); + } + + /* Need this for Lab conversions */ + icmAry2XYZ(x.twN, x.twh); + + /* Figure out the black point target */ + { + double tbL[3]; + double tbkLab[3]; + + icmXYZ2Lab(&x.twN, tbkLab, x.bk); /* Convert measured black to Lab */ + +//printf("~1 black point Lab = %f %f %f\n", tbkLab[0], tbkLab[1], tbkLab[2]); + + /* Now blend the a* b* with that of the target white point */ + tbL[0] = tbkLab[0]; + tbL[1] = bkcorrect * 0.0 + (1.0 - bkcorrect) * tbkLab[1]; + tbL[2] = bkcorrect * 0.0 + (1.0 - bkcorrect) * tbkLab[2]; + +//printf("~1 blended black Lab = %f %f %f\n", tbL[0], tbL[1], tbL[2]); + + if (bkbright > 0.0 && (bkbright <= x.bk[1] || (2.0 * bkbright) >= x.twh[1])) + warning("Black brigtness %f ignored because it is out of range",bkbright); + else if (bkbright > 0.0) { + double bkbxyz[3], bkbLab[3], vv, bl; + /* Figure out the L value of the target brightness */ + bkbxyz[0] = 0.0; + bkbxyz[1] = bkbright; + bkbxyz[2] = 0.0; + icmXYZ2Lab(&x.twN, bkbLab, bkbxyz); + + /* Do crossover to white neutral */ + vv = bkbLab[0] / (100.0 - tbL[0]); + bl = pow((1.0 - vv), x.nbrate); /* Crossover near the black */ + tbL[0] = bkbLab[0]; + tbL[1] = (1.0 - bl) * 0.0 + bl * tbL[1]; + tbL[2] = (1.0 - bl) * 0.0 + bl * tbL[2]; +//printf("~1 brighted black Lab = %f %f %f\n", tbL[0], tbL[1], tbL[2]); + } + + /* And make this the black hue to aim for */ + icmLab2XYZ(&x.twN, x.tbk, tbL); + icmAry2XYZ(x.tbN, x.tbk); + if (verb) + printf("Adjusted target black XYZ %.2f %.2f %.2f, Lab %.2f %.2f %.2f\n", + x.tbk[0], x.tbk[1], x.tbk[2], tbL[0], tbL[1], tbL[2]); + } + + /* Figure out the gamma curve black offset value */ + /* that will give us the black level we actually have. */ + { + double yy, tby; /* Target black y */ + + /* Make target black Y as high as necessary */ + /* to get the black point hue */ + /* ????? should do this by increasing L* until XYZ > x.bk ????? */ + tby = x.bk[1]; +//printf("Target Y from Y = %f\n",tby); + yy = x.bk[0] * x.tbk[1]/x.tbk[0]; +//printf("Target Y from X = %f\n",yy); + if (yy > tby) + tby = yy; + yy = x.bk[2] * x.tbk[1]/x.tbk[2]; +//printf("Target Y from Z = %f\n",yy); + if (yy > tby) + tby = yy; + + if (x.tbk[1] > tby) /* If target is already high enough */ + tby = x.tbk[1]; + + if (verb) { + double tbp[3], tbplab[3]; + tbp[0] = x.tbk[0] * tby/x.tbk[1]; + tbp[1] = tby; + tbp[2] = x.tbk[2] * tby/x.tbk[1]; + icmXYZ2Lab(&x.twN, tbplab, tbp); + printf("Target black after min adjust: XYZ %.3f %.3f %.3f, Lab %.3f %.3f %.3f\n", + tbp[0], tbp[1], tbp[2], tbplab[0], tbplab[1], tbplab[2]); + } + + /* Figure out the x.gioff and egamma needed to get this x.gooff and gamma */ + x.gooff = tby / x.twh[1]; /* Convert to relative */ + + /* tech_gamma() does the hard work */ + tech_gamma(&x, &x.egamma, &x.gooff, &x.gioff, egamma, gamma, x.gooff); + } + + if (verb) + printf("Gamma curve input offset = %f, output offset = %f, power = %f\n",x.gioff,x.gooff,x.egamma); + + /* For ambient light compensation, we make use of CIECAM02 */ + if (ambient > 0.0) { + double xyz[3], Jab[3]; + double t1, t0, a1, a0; + + /* Setup default source viewing conditions */ + if ((x.svc = new_icxcam(cam_default)) == NULL + || (x.dvc = new_icxcam(cam_default)) == NULL) { + error("Failed to create source and destination CAMs"); + } + + switch(x.gammat) { + case gt_power: /* There's nothing obvious for these cases, */ + case gt_Lab: /* So default to a computerish source viewing condition */ + + case gt_sRGB: /* sRGB standard viewing conditions */ + x.svc->set_view(x.svc, vc_none, + x.nwh, /* Display normalised white point */ + 0.2 * 80.0, /* Adapting luminence, 20% of display 80 cd/m^2 */ + 0.2, /* Background relative to reference white */ + 80.0, /* Display is 80 cd/m^2 */ + 0.01, x.nwh, /* 1% flare same white point */ + 0); + break; + + case gt_Rec709: + case gt_SMPTE240M: /* Television studio conditions */ + x.svc->set_view(x.svc, vc_none, + x.nwh, /* Display normalised white point */ + 0.2 * 1000.0/3.1415, /* Adapting luminence, 20% of 1000 lux in cd/m^2 */ + 0.2, /* Background relative to reference white */ + 1000.0/3.1415, /* Luminance of white in the Image field (cd/m^2) */ + 0.01, x.nwh, /* 1% flare same white point */ + 0); + break; + + default: + error("Unknown gamma type"); + } + /* The display we're calibratings situation */ + x.dvc->set_view(x.dvc, vc_none, + x.nwh, /* Display normalised white point */ + 0.2 * ambient, /* Adapting luminence, 20% of ambient in cd/m^2 */ + 0.2, /* Background relative to reference white */ + x.twh[1], /* Target white level (cd/m^2) */ + 0.01, x.nwh, /* 1% flare same white point */ + 0); + + /* Compute the normalisation values */ + x.svc->XYZ_to_cam(x.svc, Jab, x.nwh); /* Relative white point */ + x.dvc->cam_to_XYZ(x.dvc, xyz, Jab); + t1 = x.nwh[1]; + a1 = xyz[1]; + + xyz[0] = x.tbk[1]/x.twh[1] * x.nwh[0]; + xyz[1] = x.tbk[1]/x.twh[1] * x.nwh[1]; + xyz[2] = x.tbk[1]/x.twh[1] * x.nwh[2]; + t0 = xyz[1]; + x.svc->XYZ_to_cam(x.svc, Jab, xyz); /* Relative black Y */ + x.dvc->cam_to_XYZ(x.dvc, xyz, Jab); + a0 = xyz[1]; + +//printf("~1 t1 = %f, t0 = %f\n",t1,t0); +//printf("~1 a1 = %f, a0 = %f\n",a1,a0); + x.vn1 = (t1 - t0)/(a1 - a0); /* Scale factor */ + x.vn0 = t0 - (a0 * x.vn1); /* Then offset */ +//printf("~1 vn1 = %f, vn0 = %f\n",x.vn1, x.vn0); +//printf("~1 fix a1 = %f, should be = %f\n",a1 * x.vn1 + x.vn0, t1); +//printf("~1 fix a0 = %f, should be = %f\n",a0 * x.vn1 + x.vn0, t0); + + x.vc = 1; + + /* Compute aproximate power of viewing transform */ + if (verb) { + double v; + v = view_xform(&x, 0.5); + v = log(v) / log(0.5); + printf("Viewing conditions adjustment aprox. power = %f\n",v); + } +#ifdef NEVER +{ + int i; + + printf("~1 viewing xtranform:\n"); + for (i = 0; i <= 10; i++) { + double w, v = i/10.0; + + w = view_xform(&x, v); + printf("~1 in %f -> %f\n",v,w); + } +} +#endif /* NEVER */ + } + + /* - - - - - - - - - - - - - - - - - - - - - */ + /* Start with a scaled down number of test points and refine threshold, */ + /* and double/halve these on each iteration. */ + if (verb && verify != 2) + printf("Total Iteration %d, Final Samples = %d Final Repeat threshold = %f\n", + mxits, rsteps, errthr); + if (verify == 2) { + rsteps = VER_RES; + errthr = 0.0; + } else { + rsteps /= (1 << (mxits-1)); + errthr *= pow((double)(1 << (mxits-1)), THRESH_SCALE_POW); + } + + /* Setup the initial calibration test point values */ + init_csamp(&asgrey, &x, doupdate, verify, verify == 2 ? 1 : 0, rsteps); + + /* Calculate the initial calibration curve values */ + if (verify != 2 && !doupdate) { + int nsamp = 128; + mcvco *sdv[3]; /* Scattered data for creating mcv */ + + for (j = 0; j < 3; j++) { + if ((x.rdac[j] = new_mcv()) == NULL) { + dr->del(dr); + error("new_mcv x.rdac[%d] failed",j); + } + } + + for (j = 0; j < 3; j++) { + if ((sdv[j] = malloc(sizeof(mcvco) * rsteps)) == NULL) { + dr->del(dr); + error ("Malloc of scattered data points failed"); + } + } + + if (verb) + printf("Creating initial calibration curves...\n"); + + /* Copy the sample points */ + for (i = 0; i < rsteps; i++) { + for (j = 0; j < 3; j++) { + sdv[j][i].p = asgrey.s[i].v; + sdv[j][i].v = asgrey.s[i].rgb[j]; + sdv[j][i].w = 1.0; + } + } + if (x.nat) /* Make curve go thought white if possible */ + sdv[0][rsteps-1].w = sdv[1][rsteps-1].w = sdv[2][rsteps-1].w = 50.0; + + /* Create an initial set of RAMDAC curves */ + for (j = 0; j < 3; j++) + x.rdac[j]->fit(x.rdac[j], 0, fitord, sdv[j], rsteps, RDAC_SMOOTH); + + /* Make sure that if we are using native brightness and white point, */ + /* that the curves go to a perfect 1.0 ... */ + if (x.nat) { + for (j = 0; j < 3; j++) + x.rdac[j]->force_1(x.rdac[j], 1.0); + } + + for (j = 0; j < 3; j++) + free (sdv[j]); + } + +#ifdef DEBUG_PLOT + /* Plot the initial curves */ + if (verify != 2) { + #define XRES 255 + double xx[XRES]; + double y1[XRES]; + double y2[XRES]; + double y3[XRES]; + double rgb[3]; + for (i = 0; i < XRES; i++) { + double drgb[3], rgb[3]; + xx[i] = i/(XRES-1.0); + rgb[0] = rgb[1] = rgb[2] = xx[i]; + for (j = 0; j < 3; j++) + drgb[j] = x.rdac[j]->interp(x.rdac[j], rgb[j]); + y1[i] = drgb[0]; + y2[i] = drgb[1]; + y3[i] = drgb[2]; + } + printf("Initial ramdac curves\n"); + do_plot(xx,y1,y2,y3,XRES); + #undef XRES + } +#endif + + dr->reset_targ_w(dr); /* Reset white drift target at start of main cal. */ + + /* Now we go into the main verify & refine loop */ + for (it = verify == 2 ? mxits : 0; it < mxits || verify != 0; + rsteps *= 2, errthr /= (it < mxits) ? pow(2.0,THRESH_SCALE_POW) : 1.0, it++) { + int totmeas = 0; /* Total number of measurements in this pass */ + col set[3]; /* Variable to read one to three values from the display */ + + /* Verify pass */ + if (it >= mxits) + rsteps = VER_RES; /* Fixed verification resolution */ + else + thrfail = 0; /* Not verify pass */ + + + /* re-init asgrey if the number of test points has changed */ + reinit_csamp(&asgrey, &x, verify, (verify == 2 || it >= mxits) ? 1 : 0, rsteps); + + if (verb) { + if (it >= mxits) + printf("Doing verify pass with %d sample points\n",rsteps); + else + printf("Doing iteration %d with %d sample points and repeat threshold of %f DE\n", + it+1,rsteps, errthr); + } + /* Read and adjust each step */ + /* Do this white to black to track drift in native white point */ + for (i = rsteps-1; i >= 0; i--) { + double rpt; + double peqXYZ[3]; /* Previous steps equivalent aim point */ + double bestrgb[3]; /* In case we fail */ + double bestxyz[3]; + double prevde = 1e7; + double bestde = 1e7; + double bestdc = 1e7; + double bestpeqde = 1e7; + double besthde = 1e7; + double rgain = REFINE_GAIN; /* Scale down if lots of repeats */ + int mjac = 0; /* We measured the Jacobian */ + + /* Setup a second termination threshold chriteria based on */ + /* the delta E to the previous step point for the last pass. */ + if (i == (rsteps-1) || it < (mxits-1)) { + icmAry2Ary(peqXYZ, asgrey.s[i].tXYZ); /* Its own aim point */ + } else { + double Lab1[3], Lab2[3], Lab3[3]; + icmXYZ2Lab(&x.twN, Lab1, asgrey.s[i+1].XYZ); + icmXYZ2Lab(&x.twN, Lab2, asgrey.s[i].tXYZ); + Lab3[0] = Lab2[0]; + Lab3[1] = 0.5 * (Lab1[1] + Lab2[1]); /* Compute aim point between actual target */ + Lab3[2] = 0.5 * (Lab1[2] + Lab2[2]); /* and previous point. */ + icmLab2XYZ(&x.twN, asgrey.s[i].tXYZ, Lab3); + Lab1[0] = Lab2[0]; /* L of current target with ab of previous as 2nd threshold */ + icmLab2XYZ(&x.twN, peqXYZ, Lab1); + } + + /* Until we meet the necessary accuracy or give up */ + for (rpt = 0;rpt < MAX_RPTS; rpt++) { + int gworse = 0; /* information flag */ + double wde; /* informational */ + + set[0].r = asgrey.s[i].rgb[0]; + set[0].g = asgrey.s[i].rgb[1]; + set[0].b = asgrey.s[i].rgb[2]; + set[0].id = NULL; + + /* Read patches (no auto cr in case we repeat last patch) */ + if ((rv = dr->read(dr, set, 1, rsteps-i, rsteps, 0, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + totmeas++; + + icmAry2Ary(asgrey.s[i].pXYZ, asgrey.s[i].XYZ); /* Remember previous XYZ */ + icmAry2Ary(asgrey.s[i].XYZ, set[0].XYZ); /* Transfer current reading */ + + /* If native white and we've just measured it, */ + /* and we're not doing a verification, */ + /* adjust all the other point targets txyz to track the white. */ + if (x.nat && i == (rsteps-1) && it < mxits && asgrey.s[i].v == 1.0) { + + icmAry2Ary(x.twh, asgrey.s[i].XYZ); /* Set current white */ + icmAry2XYZ(x.twN, x.twh); /* Need this for Lab conversions */ + init_csamp_txyz(&asgrey, &x, 1); /* Recompute txyz's */ + icmAry2Ary(peqXYZ, asgrey.s[i].tXYZ); /* Fix peqXYZ */ +//printf("~1 Just reset target white to native white\n"); + if (wdrift) { /* Make sure white drift is reset on next read. */ + dr->reset_targ_w(dr); /* Reset this */ + } + } + + /* Compute the current change wanted to hit target */ + icmSub3(asgrey.s[i].deXYZ, asgrey.s[i].tXYZ, asgrey.s[i].XYZ); + asgrey.s[i].de = icmXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i].XYZ); + asgrey.s[i].peqde = icmXYZLabDE(&x.twN, peqXYZ, asgrey.s[i].XYZ); + asgrey.s[i].hde = 0.8 * asgrey.s[i].de + 0.2 * asgrey.s[i].peqde; + /* Eudclidian difference of XYZ, because this doesn't always track Lab */ + asgrey.s[i].dc = icmLabDE(asgrey.s[i].tXYZ, asgrey.s[i].XYZ); + + /* Compute change from last XYZ */ + icmSub3(asgrey.s[i].dXYZ, asgrey.s[i].XYZ, asgrey.s[i].pXYZ); + +#ifdef DEBUG + printf("\n\nTest point %d, v = %f\n",rsteps - i,asgrey.s[i].v); + printf("Current rgb %f %f %f -> XYZ %f %f %f, de %f, dc %f\n", + asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2], + asgrey.s[i].XYZ[0], asgrey.s[i].XYZ[1], asgrey.s[i].XYZ[2], + asgrey.s[i].de, asgrey.s[i].dc); + printf("Target XYZ %f %f %f, delta needed %f %f %f\n", + asgrey.s[i].tXYZ[0], asgrey.s[i].tXYZ[1], asgrey.s[i].tXYZ[2], + asgrey.s[i].deXYZ[0], asgrey.s[i].deXYZ[1], asgrey.s[i].deXYZ[2]); + if (rpt > 0) { + printf("Intended XYZ change %f %f %f, actual change %f %f %f\n", + asgrey.s[i].pdXYZ[0], asgrey.s[i].pdXYZ[1], asgrey.s[i].pdXYZ[2], + asgrey.s[i].dXYZ[0], asgrey.s[i].dXYZ[1], asgrey.s[i].dXYZ[2]); + } +#endif + + if (it < mxits) { /* Not verify, apply correction */ + int impj = 0; /* We improved the Jacobian */ + int dclip = 0; /* We clipped the new RGB */ +#ifdef ADJ_JACOBIAN + int isclipped = 0; + +#ifndef CLIP /* Check for cliping */ + /* Don't try and update the Jacobian if the */ + /* device values are going out of gamut, */ + /* and being clipped without Jac correction being aware. */ + for (j = 0; j < 3; j++) { + if (asgrey.s[i].rgb[j] <= 0.0 || asgrey.s[i].rgb[j] >= 1.0) { + isclipped = 1; + break; + } + } +#endif /* !CLIP */ + +#ifdef REMEAS_JACOBIAN + /* If the de hasn't improved, try and measure the Jacobian */ + if (it < (rsteps-1) && mjac == 0 && asgrey.s[i].de > (0.8 * prevde)) { + double dd; + if (asgrey.s[i].v < 0.5) + dd = 0.05; + else + dd= -0.05; + set[0].r = asgrey.s[i].rgb[0] + dd; + set[0].g = asgrey.s[i].rgb[1]; + set[0].b = asgrey.s[i].rgb[2]; + set[0].id = NULL; + set[1].r = asgrey.s[i].rgb[0]; + set[1].g = asgrey.s[i].rgb[1] + dd; + set[1].b = asgrey.s[i].rgb[2]; + set[1].id = NULL; + set[2].r = asgrey.s[i].rgb[0]; + set[2].g = asgrey.s[i].rgb[1]; + set[2].b = asgrey.s[i].rgb[2] + dd; + set[2].id = NULL; + + if ((rv = dr->read(dr, set, 1, rsteps-i, rsteps, 0, 0, instClamp)) != 0 + || (rv = dr->read(dr, set+1, 1, rsteps-i, rsteps, 0, 0, instClamp)) != 0 + || (rv = dr->read(dr, set+2, 1, rsteps-i, rsteps, 0, 0, instClamp)) != 0) { + dr->del(dr); + error("display read failed with '%s'\n",disprd_err(rv)); + } + totmeas += 3; + + /* Matrix organization is J[XYZ][RGB] for del RGB->del XYZ*/ + for (j = 0; j < 3; j++) { + asgrey.s[i].j[0][j] = (set[j].XYZ[0] - asgrey.s[i].XYZ[0]) / dd; + asgrey.s[i].j[1][j] = (set[j].XYZ[1] - asgrey.s[i].XYZ[1]) / dd; + asgrey.s[i].j[2][j] = (set[j].XYZ[2] - asgrey.s[i].XYZ[2]) / dd; + } + if (icmInverse3x3(asgrey.s[i].ij, asgrey.s[i].j)) { + /* Should repeat with bigger dd ? */ +//printf("~1 matrix =\n"); +//printf("~1 %f %f %f\n", asgrey.s[i].j[0][0], asgrey.s[i].j[0][1], asgrey.s[i].j[0][2]); +//printf("~1 %f %f %f\n", asgrey.s[i].j[1][0], asgrey.s[i].j[1][1], asgrey.s[i].j[1][2]); +//printf("~1 %f %f %f\n", asgrey.s[i].j[2][0], asgrey.s[i].j[2][1], asgrey.s[i].j[2][2]); + if (verb) + printf("dispcal: inverting Jacobian failed (3) - falling back\n"); + + /* Revert to the initial Jacobian */ + icmCpy3x3(asgrey.s[i].ij, asgrey.s[i].fb_ij); + } + /* Restart at the best we've had */ + if (asgrey.s[i].hde > besthde) { + asgrey.s[i].de = bestde; + asgrey.s[i].dc = bestdc; + asgrey.s[i].peqde = bestpeqde; + asgrey.s[i].hde = besthde; + asgrey.s[i].rgb[0] = bestrgb[0]; + asgrey.s[i].rgb[1] = bestrgb[1]; + asgrey.s[i].rgb[2] = bestrgb[2]; + asgrey.s[i].XYZ[0] = bestxyz[0]; + asgrey.s[i].XYZ[1] = bestxyz[1]; + asgrey.s[i].XYZ[2] = bestxyz[2]; + icmSub3(asgrey.s[i].deXYZ, asgrey.s[i].tXYZ, asgrey.s[i].XYZ); + } + mjac = 1; + impj = 1; + } +#endif /* REMEAS_JACOBIAN */ + + /* Compute a correction to the Jacobian if we can. */ + /* (Don't do this unless we have a solid previous */ + /* reading for this patch) */ + if (impj == 0 && rpt > 0 && isclipped == 0) { + double nsdrgb; /* Norm squared of pdrgb */ + double spdrgb[3]; /* Scaled previous delta rgb */ + double dXYZerr[3]; /* Error in previous prediction */ + double jadj[3][3]; /* Adjustment to Jacobian */ + double tj[3][3]; /* Temp Jacobian */ + double itj[3][3]; /* Temp inverse Jacobian */ + + /* Use Broyden's Formula */ + icmSub3(dXYZerr, asgrey.s[i].dXYZ, asgrey.s[i].pdXYZ); +//printf("~1 Jacobian error = %f %f %f\n", dXYZerr[0], dXYZerr[1], dXYZerr[2]); + nsdrgb = icmNorm3sq(asgrey.s[i].pdrgb); + /* If there was sufficient change in device values */ + /* to be above any noise: */ + if (nsdrgb >= (0.005 * 0.005)) { + icmScale3(spdrgb, asgrey.s[i].pdrgb, 1.0/nsdrgb); + icmTensMul3(jadj, dXYZerr, spdrgb); + +#ifdef DEBUG + /* Check that new Jacobian predicts previous delta XYZ */ + { + double eXYZ[3]; + + /* Make a full adjustment to temporary Jac */ + icmAdd3x3(tj, asgrey.s[i].j, jadj); + icmMulBy3x3(eXYZ, tj, asgrey.s[i].pdrgb); + icmSub3(eXYZ, eXYZ, asgrey.s[i].dXYZ); + printf("Jac check resid %f %f %f\n", eXYZ[0], eXYZ[1], eXYZ[2]); + } +#endif /* DEBUG */ + + /* Add part of our correction to actual Jacobian */ + icmScale3x3(jadj, jadj, JAC_COR_FACT); + icmAdd3x3(tj, asgrey.s[i].j, jadj); + if (icmInverse3x3(itj, tj) == 0) { /* Invert OK */ + icmCpy3x3(asgrey.s[i].j, tj); /* Use adjusted */ + icmCpy3x3(asgrey.s[i].ij, itj); + impj = 1; + } +//else printf("~1 ij failed\n"); + } +//else printf("~1 nsdrgb was below threshold\n"); + } +//else if (isclipped) printf("~1 no j update: rgb is clipped\n"); +#endif /* ADJ_JACOBIAN */ + + /* Track the best solution we've found */ + if (asgrey.s[i].hde <= besthde) { + bestde = asgrey.s[i].de; + bestdc = asgrey.s[i].dc; + bestpeqde = asgrey.s[i].peqde; + besthde = asgrey.s[i].hde; + bestrgb[0] = asgrey.s[i].rgb[0]; + bestrgb[1] = asgrey.s[i].rgb[1]; + bestrgb[2] = asgrey.s[i].rgb[2]; + bestxyz[0] = asgrey.s[i].XYZ[0]; + bestxyz[1] = asgrey.s[i].XYZ[1]; + bestxyz[2] = asgrey.s[i].XYZ[2]; + + } else if (asgrey.s[i].dc > bestdc) { + /* we got worse in Lab and XYZ ! */ + + wde = asgrey.s[i].de; + + /* If we've wandered too far, return to best we found */ + if (asgrey.s[i].hde > (3.0 * besthde)) { + asgrey.s[i].de = bestde; + asgrey.s[i].dc = bestdc; + asgrey.s[i].peqde = bestpeqde; + asgrey.s[i].hde = besthde; + asgrey.s[i].rgb[0] = bestrgb[0]; + asgrey.s[i].rgb[1] = bestrgb[1]; + asgrey.s[i].rgb[2] = bestrgb[2]; + asgrey.s[i].XYZ[0] = bestxyz[0]; + asgrey.s[i].XYZ[1] = bestxyz[1]; + asgrey.s[i].XYZ[2] = bestxyz[2]; + icmSub3(asgrey.s[i].deXYZ, asgrey.s[i].tXYZ, asgrey.s[i].XYZ); + } + + /* If the Jacobian hasn't changed, moderate the gain */ + if (impj == 0) + rgain *= 0.8; /* We might be overshooting */ + gworse = 1; + } + + /* See if we need to repeat */ + if (asgrey.s[i].de <= errthr && asgrey.s[i].peqde < errthr) { + if (verb > 1) + if (it < (mxits-1)) + printf("Point %d Delta E %f, OK\n",rsteps - i,asgrey.s[i].de); + else + printf("Point %d Delta E %f, peqDE %f, OK\n",rsteps - i,asgrey.s[i].de, asgrey.s[i].peqde); + break; /* No more retries */ + } + if ((rpt+1) >= MAX_RPTS) { + asgrey.s[i].de = bestde; /* Restore to best we found */ + asgrey.s[i].dc = bestdc; + asgrey.s[i].peqde = bestpeqde; /* Restore to best we found */ + asgrey.s[i].hde = besthde; /* Restore to best we found */ + asgrey.s[i].rgb[0] = bestrgb[0]; + asgrey.s[i].rgb[1] = bestrgb[1]; + asgrey.s[i].rgb[2] = bestrgb[2]; + asgrey.s[i].XYZ[0] = bestxyz[0]; + asgrey.s[i].XYZ[1] = bestxyz[1]; + asgrey.s[i].XYZ[2] = bestxyz[2]; + if (verb > 1) + if (it < (mxits-1)) + printf("Point %d Delta E %f, Fail\n",rsteps - i,asgrey.s[i].de); + else + printf("Point %d Delta E %f, peqDE %f, Fail\n",rsteps - i,asgrey.s[i].de,asgrey.s[i].peqde); + thrfail = 1; /* Failed to meet target */ + if (bestde > failerr) + failerr = bestde; /* Worst failed delta E */ + break; /* No more retries */ + } + if (verb > 1) { + if (gworse) + if (it < (mxits-1)) + printf("Point %d Delta E %f, Repeat (got worse)\n", rsteps - i, wde); + else + printf("Point %d Delta E %f, peqDE %f, Repeat (got worse)\n", rsteps - i, wde,asgrey.s[i].peqde); + else + if (it < (mxits-1)) + printf("Point %d Delta E %f, Repeat\n", rsteps - i,asgrey.s[i].de); + else + printf("Point %d Delta E %f, peqDE %f, Repeat\n", rsteps - i,asgrey.s[i].de,asgrey.s[i].peqde); + } + + /* Compute refinement of rgb */ + icmMulBy3x3(asgrey.s[i].pdrgb, asgrey.s[i].ij, asgrey.s[i].deXYZ); +//printf("~1 delta needed %f %f %f -> delta RGB %f %f %f\n", +//asgrey.s[i].deXYZ[0], asgrey.s[i].deXYZ[1], asgrey.s[i].deXYZ[2], +//asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]); + + /* Gain scale */ + icmScale3(asgrey.s[i].pdrgb, asgrey.s[i].pdrgb, rgain); +//printf("~1 delta RGB after gain scale %f %f %f\n", +//asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]); + +#ifdef CLIP + /* Component wise clip */ + for (j = 0; j < 3; j++) { /* Check for clip */ + if ((-asgrey.s[i].pdrgb[j]) > asgrey.s[i].rgb[j]) { + asgrey.s[i].pdrgb[j] = -asgrey.s[i].rgb[j]; + dclip = 1; + } + if (asgrey.s[i].pdrgb[j] > (1.0 - asgrey.s[i].rgb[j])) { + asgrey.s[i].pdrgb[j] = (1.0 - asgrey.s[i].rgb[j]); + dclip = 1; + } + } +#ifdef DEBUG + if (dclip) printf("delta RGB after clip %f %f %f\n", + asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]); +#endif /* DEBUG */ +#endif /* CLIP */ + /* Compute next on the basis of this one RGB */ + icmAdd3(asgrey.s[i].rgb, asgrey.s[i].rgb, asgrey.s[i].pdrgb); + + /* Save expected change in XYZ */ + icmMulBy3x3(asgrey.s[i].pdXYZ, asgrey.s[i].j, asgrey.s[i].pdrgb); +#ifdef DEBUG + printf("New rgb %f %f %f from expected del XYZ %f %f %f\n", + asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2], + asgrey.s[i].pdXYZ[0], asgrey.s[i].pdXYZ[1], asgrey.s[i].pdXYZ[2]); +#endif + } else { /* Verification, so no repeat */ + break; + } + + prevde = asgrey.s[i].de; + } /* Next repeat */ + } /* Next resolution step */ + if (verb) + printf("\n"); /* Final return for patch count */ + +#ifdef DEBUG_PLOT + /* Plot the measured response XYZ */ + { + #define XRES 256 + double xx[XRES]; + double yy[3][XRES]; + double xyz[3]; + for (i = 0; i < XRES; i++) { + xx[i] = i/(XRES-1.0); + csamp_interp(&asgrey, xyz, xx[i]); + for (j = 0; j < 3; j++) + yy[j][i] = xyz[j]; + } + printf("Measured neutral axis XYZ\n",k); + do_plot(xx,yy[0],yy[1],yy[2],XRES); + #undef XRES + } +#endif + + /* Check out the accuracy of the results: */ + { + double ctwh[3]; /* Current target white */ + icmXYZNumber ctwN; /* Same as above as XYZNumber */ + double brerr; /* Brightness error */ + double cterr; /* Color temperature delta E */ + double mnerr; /* Maximum neutral error */ + double mnv = 0.0; /* Value where maximum error is */ + double anerr; /* Average neutral error */ + double lab1[3], lab2[3]; + + /* Brightness */ + brerr = asgrey.s[asgrey.no-1].XYZ[1] - x.twh[1]; + + /* Compensate for brightness error */ + for (j = 0; j < 3; j++) + ctwh[j] = x.twh[j] * asgrey.s[asgrey.no-1].XYZ[1]/x.twh[1]; + icmAry2XYZ(ctwN, ctwh); /* Need this for Lab conversions */ + + /* Color temperature error */ + icmXYZ2Lab(&ctwN, lab1, ctwh); /* Should be 100,0,0 */ + icmXYZ2Lab(&ctwN, lab2, asgrey.s[asgrey.no-1].XYZ); + cterr = icmLabDE(lab1, lab2); + + /* check delta E of all the sample points */ + /* We're checking against our given brightness and */ + /* white point target. */ + mnerr = anerr = 0.0; + init_csamp_txyz(&asgrey, &x, 0); /* In case the targets were tweaked */ + for (i = 0; i < asgrey.no; i++) { + double err; + + /* Re-compute de in case last pass had tweaked targets */ + asgrey.s[i].de = icmXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i].XYZ); + err = asgrey.s[i].de; +//printf("RGB %.3f -> Lab %.2f %.2f %.2f, target %.2f %.2f %.2f, DE %f\n", +//asgrey.s[i].v, lab2[0], lab2[1], lab2[2], lab1[0], lab1[1], lab1[2], err); + if (err > mnerr) { + mnerr = err; + mnv = asgrey.s[i].v; + } + anerr += err; + } + anerr /= (double)asgrey.no; + + if (verb || it >= mxits) { + if (it >= mxits) + printf("Verification results:\n"); + printf("Brightness error = %f cd/m^2 (is %f, should be %f)\n",brerr,asgrey.s[asgrey.no-1].XYZ[1],x.twh[1]); + printf("White point error = %f deltaE\n",cterr); + printf("Maximum neutral error (@ %f) = %f deltaE\n",mnv, mnerr); + printf("Average neutral error = %f deltaE\n",anerr); + if (it < mxits && thrfail) + printf("Failed to meet target %f delta E, got worst case %f\n",errthr,failerr); + printf("Number of measurements taken = %d\n",totmeas); + } + } + + /* Verify loop exit */ + if (it >= (mxits + nver -1)) { + break; + } + + /* Convert our test points into calibration curves. */ + /* The call to reinit_csamp() will then convert the */ + /* curves back to current test point values. */ + /* This applies some level of cohesion between the test points, */ + /* as well as forcing monotomicity */ + if (it < mxits) { + mcvco *sdv[3]; /* Scattered data for mcv */ + + for (j = 0; j < 3; j++) { + if ((sdv[j] = malloc(sizeof(mcvco) * asgrey.no)) == NULL) { + dr->del(dr); + error ("Malloc of scattered data points failed"); + } + } + + if (verb) + printf("Computing update to calibration curves...\n"); + + /* Use fixed rgb's */ + for (j = 0; j < 3; j++) { + for (i = 0; i < asgrey.no; i++) { + sdv[j][i].p = asgrey.s[i].v; + sdv[j][i].v = asgrey.s[i].rgb[j]; + sdv[j][i].w = 1.0; +// ~~999 +#ifdef NEVER + printf("rdac %d point %d = %f, %f\n",j,i,sdv[j][i].p,sdv[j][i].v); +#endif + } + } + if (x.nat) /* Make curve go thought white if possible */ + sdv[0][rsteps-1].w = sdv[1][rsteps-1].w = sdv[2][rsteps-1].w = 10.0; + + for (j = 0; j < 3; j++) + x.rdac[j]->fit(x.rdac[j], 0, fitord, sdv[j], asgrey.no, RDAC_SMOOTH); + + /* Make sure that if we are using native brightness and white point, */ + /* that the curves go to a perfect 1.0 ... */ + if (x.nat) { + for (j = 0; j < 3; j++) + x.rdac[j]->force_1(x.rdac[j], 1.0); + } + + for (j = 0; j < 3; j++) + free(sdv[j]); +#ifdef DEBUG_PLOT + /* Plot the current curves */ + { + #define XRES 255 + double xx[XRES]; + double y1[XRES]; + double y2[XRES]; + double y3[XRES]; + double rgb[3]; + for (i = 0; i < XRES; i++) { + double drgb[3], rgb[3]; + xx[i] = i/(XRES-1.0); + rgb[0] = rgb[1] = rgb[2] = xx[i]; + for (j = 0; j < 3; j++) + drgb[j] = x.rdac[j]->interp(x.rdac[j], rgb[j]); + y1[i] = drgb[0]; + y2[i] = drgb[1]; + y3[i] = drgb[2]; + } + printf("Current ramdac curves\n"); + do_plot(xx,y1,y2,y3,XRES); + #undef XRES + } +#endif + } + } /* Next refine/verify loop */ + + free_alloc_csamp(&asgrey); /* We're done with test points */ + dr->del(dr); /* Now we're done with test window */ + + /* Write out the resulting calibration file */ + if (verify != 2) { + int calres = CAL_RES; /* steps in calibration table saved */ + cgats *ocg; /* output cgats structure */ + time_t clk = time(0); + struct tm *tsp = localtime(&clk); + char *atm = asctime(tsp); /* Ascii time */ + cgats_set_elem *setel; /* Array of set value elements */ + int ncps; /* Number of curve parameters */ + double *cps[3]; /* Arrays of curve parameters */ + char *bp = NULL, buf[100]; /* Buffer to sprintf into */ + + ocg = new_cgats(); /* Create a CGATS structure */ + ocg->add_other(ocg, "CAL"); /* our special type is Calibration file */ + + ocg->add_table(ocg, tt_other, 0); /* Add a table for RAMDAC values */ + ocg->add_kword(ocg, 0, "DESCRIPTOR", "Argyll Device Calibration Curves",NULL); + ocg->add_kword(ocg, 0, "ORIGINATOR", "Argyll dispcal", NULL); + atm[strlen(atm)-1] = '\000'; /* Remove \n from end */ + ocg->add_kword(ocg, 0, "CREATED",atm, NULL); + + ocg->add_kword(ocg, 0, "DEVICE_CLASS","DISPLAY", NULL); + ocg->add_kword(ocg, 0, "COLOR_REP","RGB", NULL); + /* Tell downstream whether they can expect that this calibration */ + /* will be applied in hardware or not. */ + ocg->add_kword(ocg, 0, "VIDEO_LUT_CALIBRATION_POSSIBLE",noramdac ? "NO" : "YES", NULL); + + /* Put the target parameters in the CGATS file too */ + if (dtype != 0) { + sprintf(buf,"%c",dtype); + ocg->add_kword(ocg, 0, "DEVICE_TYPE", buf, NULL); + } + + if (wpx == 0.0 && wpy == 0.0 && temp == 0.0 && tbright == 0.0) + ocg->add_kword(ocg, 0, "NATIVE_TARGET_WHITE","", NULL); + + sprintf(buf,"%f %f %f", x.twh[0], x.twh[1], x.twh[2]); + ocg->add_kword(ocg, 0, "TARGET_WHITE_XYZ",buf, NULL); + + switch(x.gammat) { + case gt_power: + if (egamma > 0.0) + sprintf(buf,"%f", -egamma); + else + sprintf(buf,"%f", gamma); + break; + case gt_Lab: + strcpy(buf,"L_STAR"); + break; + case gt_sRGB: + strcpy(buf,"sRGB"); + break; + case gt_Rec709: + strcpy(buf,"REC709"); + break; + case gt_SMPTE240M: + strcpy(buf,"SMPTE240M"); + break; + default: + error("Unknown gamma type"); + } + ocg->add_kword(ocg, 0, "TARGET_GAMMA",buf, NULL); + + sprintf(buf,"%f", x.oofff); + ocg->add_kword(ocg, 0, "DEGREE_OF_BLACK_OUTPUT_OFFSET",buf, NULL); + + sprintf(buf,"%f", bkcorrect); + ocg->add_kword(ocg, 0, "BLACK_POINT_CORRECTION", buf, NULL); + + sprintf(buf,"%f", x.nbrate); + ocg->add_kword(ocg, 0, "BLACK_NEUTRAL_BLEND_RATE", buf, NULL); + + if (bkbright > 0.0) { + sprintf(buf,"%f", bkbright); + ocg->add_kword(ocg, 0, "TARGET_BLACK_BRIGHTNESS",buf, NULL); + } + + /* Write rest of setup */ + switch (quality) { + case -3: /* Test value */ + bp = "ultra low"; + break; + case -2: /* Very low */ + bp = "very low"; + break; + case -1: /* Low */ + bp = "low"; + break; + case 0: /* Medum */ + bp = "medium"; + break; + case 1: /* High */ + bp = "high"; + break; + case 2: /* Ultra */ + bp = "ultra high"; + break; + default: + error("unknown quality level %d",quality); + } + ocg->add_kword(ocg, 0, "QUALITY",bp, NULL); + + ocg->add_field(ocg, 0, "RGB_I", r_t); + ocg->add_field(ocg, 0, "RGB_R", r_t); + ocg->add_field(ocg, 0, "RGB_G", r_t); + ocg->add_field(ocg, 0, "RGB_B", r_t); + + if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * 4)) == NULL) + error("Malloc failed!"); + + /* Write the video lut curve values */ + for (i = 0; i < calres; i++) { + double vv, rgb[3]; + +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(i); +#endif + vv = i/(calres-1.0); + for (j = 0; j < 3; j++) { + double cc; + cc = x.rdac[j]->interp(x.rdac[j], vv); + if (cc < 0.0) + cc = 0.0; + else if (cc > 1.0) + cc = 1.0; + rgb[j] = cc; + } + + setel[0].d = vv; + setel[1].d = rgb[0]; + setel[2].d = rgb[1]; + setel[3].d = rgb[2]; + + ocg->add_setarr(ocg, 0, setel); + } + + free(setel); + + /* Write some of the device model information to a second */ + /* table, so that we can update the calibration latter on without */ + /* having to read R,G & B curves. */ + + ocg->add_table(ocg, tt_other, 0); /* Add a second table for setup and model */ + ocg->add_kword(ocg, 1, "DESCRIPTOR", "Argyll Calibration options and model",NULL); + ocg->add_kword(ocg, 1, "ORIGINATOR", "Argyll dispcal", NULL); + atm[strlen(atm)-1] = '\000'; /* Remove \n from end */ + ocg->add_kword(ocg, 1, "CREATED",atm, NULL); + + + /* Write device model curves */ + ocg->add_field(ocg, 1, "R_P", r_t); + ocg->add_field(ocg, 1, "G_P", r_t); + ocg->add_field(ocg, 1, "B_P", r_t); + + if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * 3)) == NULL) + error("Malloc failed!"); + + ncps = -1; + for (i = 0; i < 3; i++) { + int nn; + nn = x.dcvs[i]->get_params(x.dcvs[i], &cps[i]); + if (ncps != -1 && ncps != nn) + error("Expect device model linearisation curves to have the same order"); + ncps = nn; + } + + for (i = 0; i < ncps; i++) { + setel[0].d = cps[0][i]; + setel[1].d = cps[1][i]; + setel[2].d = cps[2][i]; + ocg->add_setarr(ocg, 1, setel); + } + + for (i = 0; i < 3; i++) + free(cps[i]); + free(setel); + + if (ocg->write_name(ocg, outname)) + error("Write error : %s",ocg->err); + + if (verb) + printf("Written calibration file '%s'\n",outname); + + ocg->del(ocg); /* Clean up */ + + } + + /* Update the ICC file with the new 'vcgt' curves */ + if (verify != 2 && doupdate && doprofile) { + icmFile *ic_fp; + icc *icco; + int j, i; + icmVideoCardGamma *wo; + + if ((icco = new_icc()) == NULL) + error ("Creation of ICC object to read profile '%s' failed",iccoutname); + + /* Open up the profile for reading */ + if ((ic_fp = new_icmFileStd_name(iccoutname,"r")) == NULL) + error ("Can't open file '%s'",iccoutname); + + /* Read header etc. */ + if ((rv = icco->read(icco,ic_fp,0)) != 0) + error ("Reading profile '%s' failed with %d, %s",iccoutname, rv,icco->err); + + /* Read every tag */ + if (icco->read_all_tags(icco) != 0) { + error("Unable to read all tags from '%s': %d, %s",iccoutname, icco->errc,icco->err); + } + + ic_fp->del(ic_fp); + + wo = (icmVideoCardGamma *)icco->read_tag(icco, icSigVideoCardGammaTag); + if (wo == NULL) + error("Can't find VideoCardGamma tag in file '%s': %d, %s", + iccoutname, icco->errc,icco->err); + + wo->tagType = icmVideoCardGammaTableType; + wo->u.table.channels = 3; /* rgb */ + wo->u.table.entryCount = CAL_RES; /* full lut */ + wo->u.table.entrySize = 2; /* 16 bits */ + wo->allocate((icmBase*)wo); + for (j = 0; j < 3; j++) { + for (i = 0; i < CAL_RES; i++) { + double cc, vv; +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(i); +#endif + vv = i/(CAL_RES-1.0); + + cc = x.rdac[j]->interp(x.rdac[j], vv); + + if (cc < 0.0) + cc = 0.0; + else if (cc > 1.0) + cc = 1.0; + ((unsigned short*)wo->u.table.data)[CAL_RES * j + i] = (int)(cc * 65535.0 + 0.5); + } + } + + /* Open up the profile again writing */ + if ((ic_fp = new_icmFileStd_name(iccoutname,"w")) == NULL) + error ("Can't open file '%s' for writing",iccoutname); + + if ((rv = icco->write(icco,ic_fp,0)) != 0) + error ("Write to file '%s' failed: %d, %s",iccoutname, rv,icco->err); + + if (verb) + printf("Updated profile '%s'\n",iccoutname); + + ic_fp->del(ic_fp); + icco->del(icco); + + /* Create a fast matrix/shaper profile */ + /* + [ Another way of doing this would be to run all the + measured points through the calibration curves, and + then re-fit the curve/matrix to the calibrated points. + This might be more accurate ?] + + Ideally we should also re-measure primaries through calibration + rather than computing the calibrated values ? + + */ + + } else if (verify != 2 && doprofile) { + icmFile *wr_fp; + icc *wr_icco; + double uwp[3]; /* Absolute Uncalibrated White point in XYZ */ + double wp[3]; /* Absolute White point in XYZ */ + double bp[3]; /* Absolute Black point in XYZ */ + double mat[3][3]; /* Device to XYZ matrix */ + double calrgb[3]; /* 1.0 through calibration curves */ + double clrgb[3]; /* 1.0 through calibration and linearization */ + + /* Lookup white and black points */ + { + int j; + double rgb[3]; + + calrgb[0] = calrgb[1] = calrgb[2] = 1.0; + + fwddev(&x, uwp, calrgb); /* absolute uncalibrated WP (native white point) */ + +//printf("~1 native abs white point XYZ %f %f %f\n", uwp[0], uwp[1], uwp[2]); + + /* RGB 1.0 Through calibration */ + for (j = 0; j < 3; j++) { + calrgb[j] = x.rdac[j]->interp(x.rdac[j], calrgb[j]); + if (calrgb[j] < 0.0) + calrgb[j] = 0.0; + else if (calrgb[j] > 1.0) + calrgb[j] = 1.0; + } + fwddev(&x, wp, calrgb); /* absolute calibrated WP */ +//printf("~1 calibrated rgb = %f %f %f\n", calrgb[0], calrgb[1], calrgb[2]); +//printf("~1 calibrated abs white point XYZ %f %f %f\n", wp[0], wp[1], wp[2]); + + for (j = 0; j < 3; j++) + clrgb[j] = x.dcvs[j]->interp(x.dcvs[j], calrgb[j]); +//printf("~1 cal & lin rgb = %f %f %f\n", clrgb[0], clrgb[1], clrgb[2]); + + rgb[0] = rgb[1] = rgb[2] = 0.0; + + /* RGB 0.0 through calibration */ + for (j = 0; j < 3; j++) { + rgb[j] = x.rdac[j]->interp(x.rdac[j], rgb[j]); + if (rgb[j] < 0.0) + rgb[j] = 0.0; + else if (rgb[j] > 1.0) + rgb[j] = 1.0; + } + fwddev(&x, bp, rgb); /* Absolute calibrated BP */ + } + + /* Apply calibration to matrix, and then adjust it to be */ + /* relative to D50 white point, rather than absolute. */ + { + double rgb[3]; + icmXYZNumber swp; + + /* Transfer from parameter to matrix */ + icmCpy3x3(mat, x.fm); + + /* Compute the calibrated matrix values so that the curves */ + /* device curves end at 1.0. */ + + /* In the HW calibrated case this represents the lower XYZ due to */ + /* the HW calibrated lower RGB values of white compared to the raw */ + /* model response, so that the calibration curve concatentation with the */ + /* device curves can be scaled up to end at 1.0. */ + if (noramdac == 0) { + for (j = 0; j < 3; j++) { + for (i = 0; i < 3; i++) + rgb[i] = 0.0; + rgb[j] = clrgb[j]; + icmMulBy3x3(rgb, x.fm, rgb); /* clrgb -> matrix -> RGB */ + for (i = 0; i < 3; i++) + mat[i][j] = rgb[i]; + } +#ifdef NEVER + { + double rgb[3], xyz[3], lab[3]; + + rgb[0] = rgb[1] = rgb[2] = 1.0; + icmMulBy3x3(xyz, mat, rgb); + icmXYZ2Lab(&icmD50, lab, xyz); + + printf("RGB 1 through matrix = XYZ %f %f %f, Lab %f %f %f\n", xyz[0], xyz[1], xyz[2], lab[0], lab[1], lab[2]); + } +#endif + /* Adapt matrix */ + icmAry2XYZ(swp, wp); + icmChromAdaptMatrix(ICM_CAM_MULMATRIX | ICM_CAM_BRADFORD, icmD50, swp, mat); +#ifdef NEVER + { + double rgb[3], xyz[3], lab[3]; + + rgb[0] = rgb[1] = rgb[2] = 1.0; + icmMulBy3x3(xyz, mat, rgb); + icmXYZ2Lab(&icmD50, lab, xyz); + + printf("RGB 1 through chrom matrix = XYZ %f %f %f, Lab %f %f %f\n", xyz[0], xyz[1], xyz[2], lab[0], lab[1], lab[2]); + } +#endif + + /* For the calibration incororated in profile case, we should boost the */ + /* XYZ by 1/calrgb[] so that the lower calibrated RGB values results in the native */ + /* white point, but we want to reduce it by callinrgb[] to move from the native */ + /* white point to the calibrated white point. */ + } else { + icmCpy3x3(mat, x.fm); + + for (j = 0; j < 3; j++) { + for (i = 0; i < 3; i++) + rgb[i] = 0.0; + rgb[j] = clrgb[j]/calrgb[j]; + icmMulBy3x3(rgb, x.fm, rgb); /* 1/calrgb -> matrix -> RGB */ + for (i = 0; i < 3; i++) + mat[i][j] = rgb[i]; + } +#ifdef NEVER + { + double rgb[3], xyz[3], lab[3]; + + for (j = 0; j < 3; j++) + rgb[j] = calrgb[j]; + icmMulBy3x3(xyz, mat, rgb); + icmXYZ2Lab(&icmD50, lab, xyz); + + printf("RGB cal through matrix = XYZ %f %f %f, Lab %f %f %f\n", xyz[0], xyz[1], xyz[2], lab[0], lab[1], lab[2]); + } +#endif + /* Adapt matrix */ + icmAry2XYZ(swp, wp); + icmChromAdaptMatrix(ICM_CAM_MULMATRIX | ICM_CAM_BRADFORD, icmD50, swp, mat); +#ifdef NEVER + { + double rgb[3], xyz[3], lab[3]; + + icmMulBy3x3(xyz, mat, calrgb); + icmXYZ2Lab(&icmD50, lab, xyz); + + printf("RGB cal through chrom matrix = XYZ %f %f %f, Lab %f %f %f\n", xyz[0], xyz[1], xyz[2], lab[0], lab[1], lab[2]); + } +#endif + } + } + + /* Open up the file for writing */ + if ((wr_fp = new_icmFileStd_name(iccoutname,"w")) == NULL) + error("Write: Can't open file '%s'",iccoutname); + + if ((wr_icco = new_icc()) == NULL) + error("Write: Creation of ICC object failed"); + + /* Add all the tags required */ + + /* The header: */ + { + icmHeader *wh = wr_icco->header; + + /* Values that must be set before writing */ + wh->deviceClass = icSigDisplayClass; + wh->colorSpace = icSigRgbData; /* Display is RGB */ + wh->pcs = icSigXYZData; /* XYZ for matrix based profile */ + wh->renderingIntent = icRelativeColorimetric; /* For want of something */ + + wh->manufacturer = icmSigUnknownType; + wh->model = icmSigUnknownType; +#ifdef NT + wh->platform = icSigMicrosoft; +#endif +#ifdef __APPLE__ + wh->platform = icSigMacintosh; +#endif +#if defined(UNIX_X11) + wh->platform = icmSig_nix; +#endif + } + + /* Profile Description Tag: */ + { + icmTextDescription *wo; + char *dst, dstm[200]; /* description */ + + if (profDesc != NULL) + dst = profDesc; + else { + dst = iccoutname; + } + + if ((wo = (icmTextDescription *)wr_icco->add_tag( + wr_icco, icSigProfileDescriptionTag, icSigTextDescriptionType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = strlen(dst)+1; /* Allocated and used size of desc, inc null */ + wo->allocate((icmBase *)wo);/* Allocate space */ + strcpy(wo->desc, dst); /* Copy the string in */ + } + /* Copyright Tag: */ + { + icmText *wo; + char *crt; + + if (copyright != NULL) + crt = copyright; + else + crt = "Copyright, the creator of this profile"; + + if ((wo = (icmText *)wr_icco->add_tag( + wr_icco, icSigCopyrightTag, icSigTextType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = strlen(crt)+1; /* Allocated and used size of text, inc null */ + wo->allocate((icmBase *)wo);/* Allocate space */ + strcpy(wo->data, crt); /* Copy the text in */ + } + /* Device Manufacturers Description Tag: */ + if (deviceMfgDesc != NULL) { + icmTextDescription *wo; + char *dst = deviceMfgDesc; + + if ((wo = (icmTextDescription *)wr_icco->add_tag( + wr_icco, icSigDeviceMfgDescTag, icSigTextDescriptionType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = strlen(dst)+1; /* Allocated and used size of desc, inc null */ + wo->allocate((icmBase *)wo);/* Allocate space */ + strcpy(wo->desc, dst); /* Copy the string in */ + } + /* Model Description Tag: */ + if (modelDesc != NULL) { + icmTextDescription *wo; + char *dst = modelDesc; + + if ((wo = (icmTextDescription *)wr_icco->add_tag( + wr_icco, icSigDeviceModelDescTag, icSigTextDescriptionType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = strlen(dst)+1; /* Allocated and used size of desc, inc null */ + wo->allocate((icmBase *)wo);/* Allocate space */ + strcpy(wo->desc, dst); /* Copy the string in */ + } + /* Luminance tag */ + { + icmXYZArray *wo;; + + if ((wo = (icmXYZArray *)wr_icco->add_tag( + wr_icco, icSigLuminanceTag, icSigXYZArrayType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = 1; + wo->allocate((icmBase *)wo); /* Allocate space */ + wo->data[0].X = 0.0; + wo->data[0].Y = dispLum * wp[1]/uwp[1]; /* Adjust for effect of calibration */ + wo->data[0].Z = 0.0; + + if (verb) + printf("Luminance XYZ = %f %f %f\n", wo->data[0].X, wo->data[0].Y, wo->data[0].Z); + } + /* White Point Tag: */ + { + icmXYZArray *wo; + /* Note that tag types icSigXYZType and icSigXYZArrayType are identical */ + if ((wo = (icmXYZArray *)wr_icco->add_tag( + wr_icco, icSigMediaWhitePointTag, icSigXYZArrayType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = 1; + wo->allocate((icmBase *)wo); /* Allocate space */ + wo->data[0].X = wp[0] * 1.0/wp[1]; + wo->data[0].Y = wp[1] * 1.0/wp[1]; + wo->data[0].Z = wp[2] * 1.0/wp[1]; + + if (verb) + printf("White point XYZ = %f %f %f\n", wo->data[0].X, wo->data[0].Y, wo->data[0].Z); + } + /* Black Point Tag: */ + { + icmXYZArray *wo; + if ((wo = (icmXYZArray *)wr_icco->add_tag( + wr_icco, icSigMediaBlackPointTag, icSigXYZArrayType)) == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->size = 1; + wo->allocate((icmBase *)wo); /* Allocate space */ + wo->data[0].X = bp[0] * 1.0/wp[1]; + wo->data[0].Y = bp[1] * 1.0/wp[1]; + wo->data[0].Z = bp[2] * 1.0/wp[1]; + + if (verb) + printf("Black point XYZ = %f %f %f\n", wo->data[0].X, wo->data[0].Y, wo->data[0].Z); + } + + /* vcgt tag, if the display has an accessible VideoLUT */ + if (noramdac == 0) { + int j, i; + icmVideoCardGamma *wo; + wo = (icmVideoCardGamma *)wr_icco->add_tag(wr_icco, + icSigVideoCardGammaTag, icSigVideoCardGammaType); + if (wo == NULL) + error("add_tag failed: %d, %s",wr_icco->errc,wr_icco->err); + + wo->tagType = icmVideoCardGammaTableType; + wo->u.table.channels = 3; /* rgb */ + wo->u.table.entryCount = CAL_RES; /* full lut */ + wo->u.table.entrySize = 2; /* 16 bits */ + wo->allocate((icmBase*)wo); + for (j = 0; j < 3; j++) { + for (i = 0; i < CAL_RES; i++) { + double cc, vv = i/(CAL_RES-1.0); + cc = x.rdac[j]->interp(x.rdac[j], vv); + if (cc < 0.0) + cc = 0.0; + else if (cc > 1.0) + cc = 1.0; + ((unsigned short*)wo->u.table.data)[CAL_RES * j + i] = (int)(cc * 65535.0 + 0.5); + } + } + } + + /* Red, Green and Blue Colorant Tags: */ + { + icmXYZArray *wor, *wog, *wob; + if ((wor = (icmXYZArray *)wr_icco->add_tag( + wr_icco, icSigRedColorantTag, icSigXYZArrayType)) == NULL) + error("add_tag failed: %d, %s",rv,wr_icco->err); + if ((wog = (icmXYZArray *)wr_icco->add_tag( + wr_icco, icSigGreenColorantTag, icSigXYZArrayType)) == NULL) + error("add_tag failed: %d, %s",rv,wr_icco->err); + if ((wob = (icmXYZArray *)wr_icco->add_tag( + wr_icco, icSigBlueColorantTag, icSigXYZArrayType)) == NULL) + error("add_tag failed: %d, %s",rv,wr_icco->err); + + wor->size = wog->size = wob->size = 1; + wor->allocate((icmBase *)wor); /* Allocate space */ + wog->allocate((icmBase *)wog); + wob->allocate((icmBase *)wob); + + wor->data[0].X = mat[0][0]; wor->data[0].Y = mat[1][0]; wor->data[0].Z = mat[2][0]; + wog->data[0].X = mat[0][1]; wog->data[0].Y = mat[1][1]; wog->data[0].Z = mat[2][1]; + wob->data[0].X = mat[0][2]; wob->data[0].Y = mat[1][2]; wob->data[0].Z = mat[2][2]; + } + + /* Red, Green and Blue Tone Reproduction Curve Tags: */ + { + icmCurve *wor, *wog, *wob; + int ui; + + if ((wor = (icmCurve *)wr_icco->add_tag( + wr_icco, icSigRedTRCTag, icSigCurveType)) == NULL) + error("add_tag failed: %d, %s",rv,wr_icco->err); + if ((wog = (icmCurve *)wr_icco->add_tag( + wr_icco, icSigGreenTRCTag, icSigCurveType)) == NULL) + error("add_tag failed: %d, %s",rv,wr_icco->err); + if ((wob = (icmCurve *)wr_icco->add_tag( + wr_icco, icSigBlueTRCTag, icSigCurveType)) == NULL) + error("add_tag failed: %d, %s",rv,wr_icco->err); + + wor->flag = wog->flag = wob->flag = icmCurveSpec; + wor->size = wog->size = wob->size = 256; /* Number of entries */ + wor->allocate((icmBase *)wor); /* Allocate space */ + wog->allocate((icmBase *)wog); + wob->allocate((icmBase *)wob); + + /* For the HW calibrated case, we have lowered the matrix */ + /* values to reflect the calibrated RGB through the native */ + /* device model, so now we can scale up the comcatenation */ + /* of the calibration and linearisation curves so that */ + /* 1.0 in maps to 1.0 out. */ + if (noramdac == 0) { + + for (ui = 0; ui < wor->size; ui++) { + double in, rgb[3]; + + for (j = 0; j < 3; j++) { +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(ui); +#endif + in = (double)ui / (wor->size - 1.0); + + /* Transform through calibration curve */ + in = x.rdac[j]->interp(x.rdac[j], in); + + if (in < 0.0) + in = 0.0; + else if (in > 1.0) + in = 1.0; + + /* Trandform though device model linearisation */ + in = x.dcvs[j]->interp(x.dcvs[j], in); + + /* Scale back so that 1.0 in gets 1.0 out */ + in /= clrgb[j]; + + + if (in < 0.0) + in = 0.0; + else if (in > 1.0) + in = 1.0; + rgb[j] = in; +//printf("Step %d, Chan %d, %f -> %f\n",ui,j,(double)ui / (wor->size - 1.0),in); + } + wor->data[ui] = rgb[0]; /* Curve values 0.0 - 1.0 */ + wog->data[ui] = rgb[1]; + wob->data[ui] = rgb[2]; + } + + /* For the calibration incororated in profile case, */ + /* we bypass the inverse calibration curve if it would */ + /* result in saturation, and then scale the overall output */ + /* back by the calrgb[] value so that the overall curve */ + /* maps 1.0 to 1.0. The scaled up values in the matrix */ + /* then result in a calibrated RGB input mapping to the */ + /* PCS white point. */ + } else { + + for (ui = 0; ui < wor->size; ui++) { + double in, rgb[3]; + + for (j = 0; j < 3; j++) { +#if defined(__APPLE__) && defined(__POWERPC__) + gcc_bug_fix(ui); +#endif + in = (double)ui / (wor->size - 1.0); + + /* If within the inversion range, */ + /* transform through the inverse calibration curve. */ + if (in < calrgb[j]) { + in = x.rdac[j]->inv_interp(x.rdac[j], in); + if (in < 0.0) + in = 0.0; + else if (in > 1.0) + in = 1.0; + /* Pass through device model linearisation. */ + in = x.dcvs[j]->interp(x.dcvs[j], in); + + /* Linearly extrapolate when outside inv range */ + } else { + in /= calrgb[j]; + } + + /* Scale it back again to 0.0 to 1.0, */ + /* which is compensated for by matrix scale. */ + in *= calrgb[j]; + + if (in < 0.0) + in = 0.0; + else if (in > 1.0) + in = 1.0; + rgb[j] = in; +//printf("Step %d, Chan %d, %f -> %f\n",ui,j,(double)ui / (wor->size - 1.0),in); + } + wor->data[ui] = rgb[0]; /* Curve values 0.0 - 1.0 */ + wog->data[ui] = rgb[1]; + wob->data[ui] = rgb[2]; + } + } + } + + /* Write the file (including all tags) out */ + if ((rv = wr_icco->write(wr_icco,wr_fp,0)) != 0) { + error("Write file: %d, %s",rv,wr_icco->err); + } + + if (verb) + printf("Created fast shaper/matrix profile '%s'\n",iccoutname); + + /* Close the file */ + wr_icco->del(wr_icco); + wr_fp->del(wr_fp); + } + + if (verify != 2) { + for (j = 0; j < 3; j++) + x.rdac[j]->del(x.rdac[j]); + + for (k = 0; k < 3; k++) + x.dcvs[k]->del(x.dcvs[k]); + } + + if (x.svc != NULL) { + x.svc->del(x.svc); + x.svc = NULL; + } + if (x.dvc != NULL) { + x.dvc->del(x.svc); + x.dvc = NULL; + } + + free_a_disppath(disp); + + return 0; +} + + -- cgit v1.2.3