--- ray/src/px/macbethcal.c 1995/10/11 11:57:41 2.2 +++ ray/src/px/macbethcal.c 1995/10/17 10:24:21 2.5 @@ -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,12 +75,21 @@ 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}; -#define NMBMOD 3 /* Number of MacBeth moderate colors */ -short mbmod[NMBMOD] = {1,2,21}; -#define NMBSAT 6 /* Number of MacBeth saturated colors */ -short mbsat[NMBSAT] = {14,12,13,15,16,17}; +short mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White}; +#define NMBMOD 17 /* Number of MacBeth unsaturated colors */ +short mbmod[NMBMOD] = { + DarkSkin,LightSkin,BlueSky,Foliage,BlueFlower,BluishGreen, + Orange,PurplishBlue,ModerateRed,Purple,YellowGreen, + Black,Neutral35,Neutral5,Neutral65,Neutral8,White + }; + +#define NMBSAT 8 /* Number of MacBeth saturated colors */ +short mbsat[NMBSAT] = { + Red,Green,Blue,Magenta,Yellow,Cyan, + Orange,OrangeYellow + }; + int xmax, ymax; /* input image dimensions */ int bounds[4][2]; /* image coordinates of chart corners */ double imgxfm[3][3]; /* coordinate transformation matrix */ @@ -65,7 +99,7 @@ COLOR picRGB[24]; /* picture colors */ COLOR bramp[NMBNEU][2]; /* brightness ramp (per primary) */ double solmat[3][3]; /* color mapping matrix */ -FILE *debugfp = NULL; +FILE *debugfp = NULL; /* debug output picture */ char *progname; extern char *malloc(); @@ -126,6 +160,9 @@ char **argv; init(); /* initialize */ getcolors(); /* get picture colors */ compute(); /* compute color mapping */ + /* print comment */ + printf("{ Color correction file computed by %s }\n", progname); + printf("{ from scanned MacBetch color chart %s }\n", argv[1]); putmapping(); /* put out color mapping */ putdebug(); /* put out debug picture */ exit(0); @@ -231,18 +268,21 @@ getcolors() /* load in picture colors */ bresp(y, x) /* piecewise linear interpolation of primaries */ COLOR y, x; { + double cv[3]; register int i, n; for (i = 0; i < 3; i++) { - n = NMBNEU; - while (n > 0 && colval(x,i) < colval(bramp[--n][0],i)) - ; - colval(y,i) = ((colval(bramp[n+1][0],i) - colval(x,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)) * colval(bramp[n+1][1],i)) / (colval(bramp[n+1][0],i) - colval(bramp[n][0],i)); + if (cv[i] < 0.) cv[i] = 0.; } + setcolor(y, cv[0], cv[1], cv[2]); } @@ -304,15 +344,15 @@ putmapping() /* put out color mapping for pcomb -f * } /* print color mapping */ printf("ro = %g*rn + %g*gn + %g*bn ;\n", - solmat[0][0], solmat[1][0], solmat[2][0]); + solmat[0][0], solmat[0][1], solmat[0][2]); printf("go = %g*rn + %g*gn + %g*bn ;\n", - solmat[0][1], solmat[1][1], solmat[2][1]); + solmat[1][0], solmat[1][1], solmat[1][2]); printf("bo = %g*rn + %g*gn + %g*bn ;\n", - solmat[0][2], solmat[1][2], solmat[2][2]); + solmat[2][0], solmat[2][1], solmat[2][2]); } -compsoln(cin, cout, n) /* solve 3x3 system */ +compsoln(cin, cout, n) /* solve 3xN system using least-squares */ COLOR cin[], cout[]; int n; { @@ -320,15 +360,28 @@ int n; double mat[3][3], invmat[3][3]; double det; double colv[3], rowv[3]; - register int i, j; + register int i, j, k; - if (n != 3) { + if (n < 3 | n > NMBMOD) { fprintf(stderr, "%s: inconsistent code!\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", @@ -342,11 +395,19 @@ int n; for (j = 0; j < 3; j++) invmat[i][j] /= det; for (i = 0; i < 3; 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++) - rowv[j] = colval(cout[j],i); - mx3d_transform(rowv, invmat, colv); - for (j = 0; j < 3; j++) - solmat[j][i] = colv[j]; + solmat[i][j] = rowv[j]; } } @@ -357,13 +418,13 @@ COLOR cout, cin; double r, g, b; bresp(cout, cin); - r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[1][0] - + colval(cout,2)*solmat[2][0]; + r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[0][1] + + colval(cout,2)*solmat[0][2]; if (r < 0) r = 0; - g = colval(cout,0)*solmat[0][1] + colval(cout,1)*solmat[1][1] - + colval(cout,2)*solmat[2][1]; + g = colval(cout,0)*solmat[1][0] + colval(cout,1)*solmat[1][1] + + colval(cout,2)*solmat[1][2]; if (g < 0) g = 0; - b = colval(cout,0)*solmat[0][2] + colval(cout,1)*solmat[1][2] + b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1] + colval(cout,2)*solmat[2][2]; if (b < 0) b = 0; setcolor(cout, r, g, b); @@ -391,6 +452,7 @@ putdebug() /* put out debugging picture */ exit(1); } /* finish debug header */ + fputformat(COLRFMT, debugfp); putc('\n', debugfp); fprtresolu(xmax, ymax, debugfp); for (y = ymax-1; y >= 0; y--) {