--- ray/src/px/macbethcal.c 1995/10/11 18:52:00 2.3 +++ ray/src/px/macbethcal.c 1995/10/20 15:26:11 2.7 @@ -18,8 +18,33 @@ static char SCCSid[] = "$SunId$ LBL"; #include "resolu.h" #include "pmap.h" - /* MacBeth colors (CIE 1931, absolute white) */ - /* computed from spectral measurements */ + /* MacBeth colors */ +#define DarkSkin 0 +#define LightSkin 1 +#define BlueSky 2 +#define Foliage 3 +#define BlueFlower 4 +#define BluishGreen 5 +#define Orange 6 +#define PurplishBlue 7 +#define ModerateRed 8 +#define Purple 9 +#define YellowGreen 10 +#define OrangeYellow 11 +#define Blue 12 +#define Green 13 +#define Red 14 +#define Yellow 15 +#define Magenta 16 +#define Cyan 17 +#define White 18 +#define Neutral8 19 +#define Neutral65 20 +#define Neutral5 21 +#define Neutral35 22 +#define Black 23 + /* computed from 5nm spectral measurements */ + /* CIE 1931 2 degree obs, equal-energy white */ float mbxyY[24][3] = { {0.462, 0.3769, 0.0932961}, /* DarkSkin */ {0.4108, 0.3542, 0.410348}, /* LightSkin */ @@ -50,22 +75,38 @@ float mbxyY[24][3] = { COLOR mbRGB[24]; /* MacBeth RGB values */ #define NMBNEU 6 /* Number of MacBeth neutral colors */ -short mbneu[NMBNEU] = {23,22,21,20,19,18}; +short mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White}; + #define NMBMOD 16 /* Number of MacBeth unsaturated colors */ -short mbmod[NMBMOD] = {0,1,2,3,4,5,6,7,8,9,10,11,19,20,21,22}; -#define NMBSAT 6 /* Number of MacBeth saturated colors */ -short mbsat[NMBSAT] = {14,12,13,15,16,17}; +short mbmod[NMBMOD] = { + DarkSkin,LightSkin,BlueSky,Foliage,BlueFlower,BluishGreen, + PurplishBlue,ModerateRed,YellowGreen,OrangeYellow, + Black,Neutral35,Neutral5,Neutral65,Neutral8,White + }; +#define NMBSAT 8 /* Number of MacBeth saturated colors */ +short mbsat[NMBSAT] = { + Red,Green,Blue,Magenta,Yellow,Cyan, + Orange,Purple + }; + +#define REQFLGS (1L< 2 && !strcmp(argv[1], "-d")) { /* debug output */ - if ((debugfp = fopen(argv[2], "w")) == NULL) { - perror(argv[2]); - exit(1); - } + for (i = 1; i < argc && argv[i][0] == '-'; i++) + switch (argv[i][1]) { + case 'd': /* debug output */ + i++; + if (badarg(argc-i, argv+i, "s")) + goto userr; + if ((debugfp = fopen(argv[i], "w")) == NULL) { + perror(argv[i]); + exit(1); + } #ifdef MSDOS - setmode(fileno(debugfp), O_BINARY); + setmode(fileno(debugfp), O_BINARY); #endif - newheader("RADIANCE", debugfp); - printargs(argc, argv, debugfp); - argv += 2; - argc -= 2; - } - if (argc != 3 && argc != 11) - goto userr; - if (strcmp(argv[1], "-") && freopen(argv[1], "r", stdin) == NULL) { + newheader("RADIANCE", debugfp); /* start */ + printargs(argc, argv, debugfp); /* header */ + break; + case 'p': /* picture position */ + if (badarg(argc-i-1, argv+i+1, "iiiiiiii")) + goto userr; + bounds[0][0] = atoi(argv[++i]); + bounds[0][1] = atoi(argv[++i]); + bounds[1][0] = atoi(argv[++i]); + bounds[1][1] = atoi(argv[++i]); + bounds[2][0] = atoi(argv[++i]); + bounds[2][1] = atoi(argv[++i]); + bounds[3][0] = atoi(argv[++i]); + bounds[3][1] = atoi(argv[++i]); + inpispic = 2; + break; + case 'c': /* color input */ + inpispic = 0; + break; + default: + goto userr; + } + /* open files */ + if (i < argc && freopen(argv[i], "r", stdin) == NULL) { perror(argv[1]); exit(1); } - if (strcmp(argv[2], "-") && freopen(argv[2], "w", stdout) == NULL) { + if (i+1 < argc && freopen(argv[i+1], "w", stdout) == NULL) { perror(argv[2]); exit(1); } + if (inpispic) { /* load input picture header */ #ifdef MSDOS - setmode(fileno(stdin), O_BINARY); + setmode(fileno(stdin), O_BINARY); #endif - if (checkheader(stdin, COLRFMT, NULL) < 0 || - fgetresolu(&xmax, &ymax, stdin) < 0) { - fprintf(stderr, "%s: bad input picture\n", progname); - exit(1); - } - /* get chart boundaries */ - if (argc == 11) { - for (i = 0; i < 4; i++) { - if (!isint(argv[2*i+3]) | !isint(argv[2*i+4])) - goto userr; - bounds[i][0] = atoi(argv[2*i+3]); - bounds[i][1] = atoi(argv[2*i+4]); + if (checkheader(stdin, COLRFMT, NULL) < 0 || + fgetresolu(&xmax, &ymax, stdin) < 0) { + fprintf(stderr, "%s: bad input picture\n", progname); + exit(1); } - } else { + } else { /* else set default xmax and ymax */ + xmax = 512; + ymax = 2*512/3; + } + if (inpispic != 2) { /* use default boundaries */ bounds[0][0] = bounds[2][0] = .029*xmax + .5; bounds[0][1] = bounds[1][1] = .956*ymax + .5; bounds[1][0] = bounds[3][0] = .971*xmax + .5; bounds[2][1] = bounds[3][1] = .056*ymax + .5; } init(); /* initialize */ - getcolors(); /* get picture colors */ + if (inpispic) /* get picture colors */ + getpicture(); + else + getcolors(); compute(); /* compute color mapping */ + /* print comment */ + printf("{ Color correction file computed by:\n\t"); + printargs(argc, argv, stdout); + printf("}\n"); putmapping(); /* put out color mapping */ - putdebug(); /* put out debug picture */ + if (debugfp != NULL) /* put out debug picture */ + if (inpispic) + picdebug(); + else + clrdebug(); exit(0); userr: - fprintf(stderr, "Usage: %s [-d dbg.pic] input.pic output.cal [xul yul xur yur xll yll xlr ylr]\n", + fprintf(stderr, +"Usage: %s [-d dbg.pic][-p xul yul xur yur xll yll xlr ylr] input.pic [output.cal]\n", progname); + fprintf(stderr, " or: %s [-d dbg.pic] -c [xyY.dat [output.cal]]\n", + progname); exit(1); } @@ -139,6 +213,7 @@ userr: init() /* initialize */ { double quad[4][2]; + register int i; /* make coordinate transformation */ quad[0][0] = bounds[0][0]; quad[0][1] = bounds[0][1]; @@ -153,16 +228,21 @@ init() /* initialize */ fprintf(stderr, "%s: bad chart boundaries\n", progname); exit(1); } + /* map MacBeth colors to RGB space */ + for (i = 0; i < 24; i++) + xyY2RGB(mbRGB[i], mbxyY[i]); } int -chartndx(x, y) /* find color number for position */ +chartndx(x, y, cvg) /* find color number for position */ int x, y; +double cvg; { double ipos[3], cpos[3]; int ix, iy; double fx, fy; + double cmin, cmax; ipos[0] = x; ipos[1] = y; @@ -176,13 +256,15 @@ int x, y; iy = cpos[1]; fx = cpos[0] - ix; fy = cpos[1] - iy; - if (fx < .35 || fx >= .65 || fy < .35 || fy >= .65) + cmin = .5*(1.-cvg); + cmax = 1. - cmin; + if (fx < cmin || fx >= cmax || fy < cmin || fy >= cmax) return(-1); return(iy*6 + ix); } -getcolors() /* load in picture colors */ +getpicture() /* load in picture colors */ { COLR *scanln; COLOR pval; @@ -197,7 +279,7 @@ getcolors() /* load in picture colors */ exit(1); } for (i = 0; i < 24; i++) { - setcolor(picRGB[i], 0., 0., 0.); + setcolor(inpRGB[i], 0., 0., 0.); ccount[i] = 0; } for (y = ymax-1; y >= 0; y--) { @@ -207,24 +289,70 @@ getcolors() /* load in picture colors */ exit(1); } for (x = 0; x < xmax; x++) { - i = chartndx(x, y); + i = chartndx(x, y, CENTCVG); if (i >= 0) { colr_color(pval, scanln[x]); - addcolor(picRGB[i], pval); + addcolor(inpRGB[i], pval); ccount[i]++; } } } - for (i = 0; i < 24; i++) { - if (ccount[i] == 0) { - fprintf(stderr, "%s: bad chart boundaries\n", + for (i = 0; i < 24; i++) { /* compute averages */ + if (ccount[i] == 0) + continue; + d = 1./ccount[i]; + scalecolor(inpRGB[i], d); + inpflags |= 1L< 24 || + fgetval(stdin, 'f', &xyYin[0]) != 1 || + fgetval(stdin, 'f', &xyYin[1]) != 1 || + fgetval(stdin, 'f', &xyYin[2]) != 1 || + xyYin[0] < 0. | xyYin[0] > 1. | + xyYin[1] < 0. | xyYin[1] > 1.) { + fprintf(stderr, "%s: bad color input data\n", progname); exit(1); } - d = 1.0/ccount[i]; - scalecolor(picRGB[i], d); + if (n == 0) { /* calibration white */ + xyY2RGB(whiteclr, xyYin); + gotwhite++; + } else { /* standard color */ + n--; + xyY2RGB(inpRGB[n], xyYin); + inpflags |= 1L< 0 && colval(x,i) < colval(bramp[--n][0],i)) - ; + for (n = 0; n < NMBNEU-2; n++) + if (colval(x,i) < colval(bramp[n+1][0],i)) + break; cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) * colval(bramp[n][1],i) + (colval(x,i) - colval(bramp[n][0],i)) * @@ -252,28 +380,26 @@ COLOR y, x; compute() /* compute color mapping */ { COLOR clrin[NMBMOD], clrout[NMBMOD]; - COLOR ctmp; - double d; - register int i; - /* map MacBeth colors to RGB space */ - for (i = 0; i < 24; i++) { - d = mbxyY[i][2] / mbxyY[i][1]; - ctmp[0] = mbxyY[i][0] * d; - ctmp[1] = mbxyY[i][2]; - ctmp[2] = (1. - mbxyY[i][0] - mbxyY[i][1]) * d; - cie_rgb(mbRGB[i], ctmp); + register int i, n; + /* did we get what we need? */ + if ((inpflags & REQFLGS) != REQFLGS) { + fprintf(stderr, "%s: missing required input colors\n", + progname); + exit(1); } /* compute piecewise luminance curve */ for (i = 0; i < NMBNEU; i++) { - copycolor(bramp[i][0], picRGB[mbneu[i]]); + copycolor(bramp[i][0], inpRGB[mbneu[i]]); copycolor(bramp[i][1], mbRGB[mbneu[i]]); } /* compute color matrix */ - for (i = 0; i < NMBMOD; i++) { - bresp(clrin[i], picRGB[mbmod[i]]); - copycolor(clrout[i], mbRGB[mbmod[i]]); - } - compsoln(clrin, clrout, NMBMOD); + for (n = 0, i = 0; i < NMBMOD; i++) + if (inpflags & 1L< NMBMOD) { + fprintf(stderr, "%s: bad number of colors to match\n", progname); exit(1); } - for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - mat[i][j] = colval(cin[j],i); + if (n == 3) + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + mat[i][j] = colval(cin[j],i); + else { /* compute A^t A */ + for (i = 0; i < 3; i++) + for (j = i; j < 3; j++) { + mat[i][j] = 0.; + for (k = 0; k < n; k++) + mat[i][j] += colval(cin[k],i) * + colval(cin[k],j); + } + for (i = 1; i < 3; i++) /* using symmetry */ + for (j = 0; j < i; j++) + mat[i][j] = mat[j][i]; + } det = mx3d_adjoint(mat, invmat); if (fabs(det) < 1e-4) { fprintf(stderr, "%s: cannot compute color mapping\n", @@ -346,60 +484,23 @@ int n; for (j = 0; j < 3; j++) invmat[i][j] /= det; for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) - colv[j] = colval(cout[j],i); + if (n == 3) + for (j = 0; j < 3; j++) + colv[j] = colval(cout[j],i); + else + for (j = 0; j < 3; j++) { + colv[j] = 0.; + for (k = 0; k < n; k++) + colv[j] += colval(cout[k],i) * + colval(cin[k],j); + } mx3d_transform(colv, invmat, rowv); for (j = 0; j < 3; j++) solmat[i][j] = rowv[j]; } } -#else -compsoln(cin, cout, n) /* solve 3xN system (N > 3) */ -COLOR cin[], cout[]; -int n; -{ - double *au[NMBMOD], *v[3], vv[3][3], auv[NMBMOD][3], w[3]; - double b[NMBMOD]; - register int i, j; - if (n > NMBMOD) { - fprintf(stderr, "%s: inconsistent code!\n", progname); - exit(1); - } - for (i = 0; i < n; i++) /* assign rectangular matrix A */ - for (j = 0; j < 3; j++) - auv[i][j] = colval(cin[i],j); - /* svdcmp indexing requires pointer offsets */ - for (j = 0; j < 3; j++) - v[j] = vv[j] - 1; - for (i = 0; i < n; i++) - au[i] = auv[i] - 1; - /* compute singular value decomposition */ -fprintf(stderr, "A:\n"); -for (i = 1; i <= n; i++) -fprintf(stderr, "%g %g %g\n", (au-1)[i][1], (au-1)[i][2], (au-1)[i][3]); - svdcmp(au-1, n, 3, w-1, v-1); -fprintf(stderr, "U:\n"); -for (i = 0; i < n; i++) -fprintf(stderr, "%g %g %g\n", auv[i][0], auv[i][1], auv[i][2]); -fprintf(stderr, "V:\n"); -for (i = 0; i < 3; i++) -fprintf(stderr, "%g %g %g\n", vv[i][0], vv[i][1], vv[i][2]); -fprintf(stderr, "W: %g %g %g\n", w[0], w[1], w[2]); - /* zero out small weights */ - for (j = 0; j < 3; j++) - if (w[j] < 1e-4) - w[j] = 0.; - /* back substitution for each row vector */ - for (j = 0; j < 3; j++) { - for (i = 0; i < n; i++) - b[i] = colval(cout[i],j); - svbksb(au-1, w-1, v-1, n, 3, b-1, solmat[j]-1); - } -} -#endif - cvtcolor(cout, cin) /* convert color according to our mapping */ COLOR cout, cin; { @@ -419,14 +520,27 @@ COLOR cout, cin; } -putdebug() /* put out debugging picture */ +xyY2RGB(rgbout, xyYin) /* convert xyY to RGB */ +COLOR rgbout; +register float xyYin[3]; { + COLOR ctmp; + double d; + + d = xyYin[2] / xyYin[1]; + ctmp[0] = xyYin[0] * d; + ctmp[1] = xyYin[2]; + ctmp[2] = (1. - xyYin[0] - xyYin[1]) * d; + cie_rgb(rgbout, ctmp); +} + + +picdebug() /* put out debugging picture */ +{ COLOR *scan; int y; register int x, i; - if (debugfp == NULL) - return; if (fseek(stdin, 0L, 0) == EOF) { fprintf(stderr, "%s: cannot seek on input picture\n", progname); exit(1); @@ -440,8 +554,10 @@ putdebug() /* put out debugging picture */ exit(1); } /* finish debug header */ + fputformat(COLRFMT, debugfp); putc('\n', debugfp); fprtresolu(xmax, ymax, debugfp); + /* write debug picture */ for (y = ymax-1; y >= 0; y--) { if (freadscan(scan, xmax, stdin) < 0) { fprintf(stderr, "%s: error rereading input picture\n", @@ -449,7 +565,7 @@ putdebug() /* put out debugging picture */ exit(1); } for (x = 0; x < xmax; x++) { - i = chartndx(x, y); + i = chartndx(x, y, CENTCVG); if (i < 0) cvtcolor(scan[x], scan[x]); else @@ -461,5 +577,57 @@ putdebug() /* put out debugging picture */ exit(1); } } + /* clean up */ + fclose(debugfp); + free((char *)scan); +} + + +clrdebug() /* put out debug picture from color input */ +{ + static COLR blkclr = BLKCOLR; + COLR mbclr[24], cvclr[24]; + COLR *scan; + COLOR ctmp; + int y; + register int i, x; + /* convert colors */ + for (i = 0; i < 24; i++) { + setcolr(mbclr[i], colval(mbRGB[i],RED), + colval(mbRGB[i],GRN), colval(mbRGB[i],BLU)); + if (inpflags & 1L<= 0; y--) { + for (x = 0; x < xmax; x++) + if ((i = chartndx(x, y, CENTCVG)) >= 0) + copycolr(scan[x], mbclr[i]); + else if ((i = chartndx(x, y, FULLCVG)) >= 0 && + inpflags & 1L<