/* * 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; }