--- ray/src/px/macbethcal.c	1995/10/20 15:26:11	2.7
+++ ray/src/px/macbethcal.c	1995/11/01 09:04:16	2.10
@@ -77,31 +77,31 @@ COLOR	mbRGB[24];		/* MacBeth RGB values */
 #define	NMBNEU		6	/* Number of MacBeth neutral colors */
 short	mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White};
 
-#define NMBMOD		16	/* Number of MacBeth unsaturated colors */
-short	mbmod[NMBMOD] = {
-		DarkSkin,LightSkin,BlueSky,Foliage,BlueFlower,BluishGreen,
-		PurplishBlue,ModerateRed,YellowGreen,OrangeYellow,
-		Black,Neutral35,Neutral5,Neutral65,Neutral8,White
-	};
+#define  NEUFLGS	(1L<<White|1L<<Neutral8|1L<<Neutral65| \
+				1L<<Neutral5|1L<<Neutral35|1L<<Black)
 
-#define NMBSAT		8	/* Number of MacBeth saturated colors */
-short	mbsat[NMBSAT] = {
-		Red,Green,Blue,Magenta,Yellow,Cyan,
-		Orange,Purple
-	};
+#define  SATFLGS	(1L<<Red|1L<<Green|1L<<Blue|1L<<Magenta|1L<<Yellow| \
+			1L<<Cyan|1L<<Orange|1L<<Purple|1L<<PurplishBlue| \
+			1L<<YellowGreen|1<<OrangeYellow|1L<<BlueFlower)
 
-#define  REQFLGS	(1L<<White|1L<<Neutral8|1L<<Neutral65| \
-				1L<<Neutral5|1L<<Neutral35|1L<<Black)
+#define  UNSFLGS	(1L<<DarkSkin|1L<<LightSkin|1L<<BlueSky|1L<<Foliage| \
+			1L<<BluishGreen|1L<<ModerateRed)
 
-#define  CENTCVG	0.3	/* measured coverage of square sample */
-#define  FULLCVG	0.9	/* coverage of entire square */
+#define  REQFLGS	NEUFLGS			/* need these colors */
+#define  MODFLGS	(NEUFLGS|UNSFLGS)	/* should be in gamut */
 
+#define  RG_BORD	0	/* patch border */
+#define  RG_CENT	01	/* central region of patch */
+#define  RG_ORIG	02	/* original color region */
+#define  RG_CORR	04	/* corrected color region */
+
 int	xmax, ymax;		/* input image dimensions */
 int	bounds[4][2];		/* image coordinates of chart corners */
 double	imgxfm[3][3];		/* coordinate transformation matrix */
 
 COLOR	inpRGB[24];		/* measured or scanned input colors */
 long	inpflags = 0;		/* flags of which colors were input */
+long	gmtflags = 0;		/* flags of out-of-gamut colors */
 
 COLOR	bramp[NMBNEU][2];	/* brightness ramp (per primary) */
 double	solmat[3][3];		/* color mapping matrix */
@@ -190,8 +190,10 @@ char	**argv;
 		getcolors();
 	compute();			/* compute color mapping */
 					/* print comment */
-	printf("{ Color correction file computed by:\n\t");
+	printf("{\n\tColor correction file computed by:\n\t\t");
 	printargs(argc, argv, stdout);
+	printf("\n\tUsage: pcomb -f %s uncorrected.pic > corrected.pic\n",
+			i+1 < argc ? argv[i+1] : "{this_file}");
 	printf("}\n");
 	putmapping();			/* put out color mapping */
 	if (debugfp != NULL)		/* put out debug picture */
@@ -235,14 +237,13 @@ init()				/* initialize */
 
 
 int
-chartndx(x, y, cvg)			/* find color number for position */
+chartndx(x, y, np)			/* find color number for position */
 int	x, y;
-double	cvg;
+int	*np;
 {
 	double	ipos[3], cpos[3];
 	int	ix, iy;
 	double	fx, fy;
-	double  cmin, cmax;
 
 	ipos[0] = x;
 	ipos[1] = y;
@@ -251,16 +252,19 @@ double	cvg;
 	cpos[0] /= cpos[2];
 	cpos[1] /= cpos[2];
 	if (cpos[0] < 0. || cpos[0] >= 6. || cpos[1] < 0. || cpos[1] >= 4.)
-		return(-1);
+		return(RG_BORD);
 	ix = cpos[0];
 	iy = cpos[1];
 	fx = cpos[0] - ix;
 	fy = cpos[1] - iy;
-	cmin = .5*(1.-cvg);
-	cmax = 1. - cmin;
-	if (fx < cmin || fx >= cmax || fy < cmin || fy >= cmax)
-		return(-1);
-	return(iy*6 + ix);
+	*np = iy*6 + ix;
+	if (fx >= 0.35 && fx < 0.65 && fy >= 0.35 && fy < 0.65)
+		return(RG_CENT);
+	if (fx < 0.05 || fx >= 0.95 || fy < 0.05 || fy >= 0.95)
+		return(RG_BORD);
+	if (fx >= 0.5)			/* right side is corrected */
+		return(RG_CORR);
+	return(RG_ORIG);		/* left side is original */
 }
 
 
@@ -270,8 +274,8 @@ getpicture()				/* load in picture colors */
 	COLOR	pval;
 	int	ccount[24];
 	double	d;
-	int	y;
-	register int	x, i;
+	int	y, i;
+	register int	x;
 
 	scanln = (COLR *)malloc(xmax*sizeof(COLR));
 	if (scanln == NULL) {
@@ -288,14 +292,12 @@ getpicture()				/* load in picture colors */
 					progname);
 			exit(1);
 		}
-		for (x = 0; x < xmax; x++) {
-			i = chartndx(x, y, CENTCVG);
-			if (i >= 0) {
+		for (x = 0; x < xmax; x++)
+			if (chartndx(x, y, &i) == RG_CENT) {
 				colr_color(pval, scanln[x]);
 				addcolor(inpRGB[i], pval);
 				ccount[i]++;
 			}
-		}
 	}
 	for (i = 0; i < 24; i++) {		/* compute averages */
 		if (ccount[i] == 0)
@@ -359,28 +361,27 @@ getcolors()			/* get xyY colors from standard input */
 bresp(y, x)		/* piecewise linear interpolation of primaries */
 COLOR	y, x;
 {
-	double	cv[3];
 	register int	i, n;
 
 	for (i = 0; i < 3; 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(y,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]);
 }
 
 
 compute()			/* compute color mapping */
 {
-	COLOR	clrin[NMBMOD], clrout[NMBMOD];
-	register int	i, n;
+	COLOR	clrin[24], clrout[24];
+	long	cflags;
+	COLOR	ctmp;
+	register int	i, j, n;
 					/* did we get what we need? */
 	if ((inpflags & REQFLGS) != REQFLGS) {
 		fprintf(stderr, "%s: missing required input colors\n",
@@ -392,14 +393,33 @@ compute()			/* compute color mapping */
 		copycolor(bramp[i][0], inpRGB[mbneu[i]]);
 		copycolor(bramp[i][1], mbRGB[mbneu[i]]);
 	}
-					/* compute color matrix */
-	for (n = 0, i = 0; i < NMBMOD; i++)
-		if (inpflags & 1L<<mbmod[i]) {
-			bresp(clrin[n], inpRGB[mbmod[i]]);
-			copycolor(clrout[n], mbRGB[mbmod[i]]);
-			n++;
-		}
-	compsoln(clrin, clrout, n);
+					/* compute color mapping */
+	do {
+		cflags = inpflags & ~gmtflags;
+		n = 0;				/* compute transform matrix */
+		for (i = 0; i < 24; i++)
+			if (cflags & 1L<<i) {
+				bresp(clrin[n], inpRGB[i]);
+				copycolor(clrout[n], mbRGB[i]);
+				n++;
+			}
+		compsoln(clrin, clrout, n);
+						/* check out-of-gamut colors */
+		for (i = 0; i < 24; i++)
+			if (cflags & 1L<<i) {
+				cresp(ctmp, mbRGB[i]);
+				for (j = 0; j < 3; j++)
+					if (colval(ctmp,j) <= 0. ||
+						colval(ctmp,j) >= 1.) {
+						gmtflags |= 1L<<i;
+						break;
+					}
+			}
+	} while (cflags & gmtflags);
+	if (gmtflags & MODFLGS)
+		fprintf(stderr,
+		"%s: warning - some moderate colors are out of gamut\n",
+				progname);
 }
 
 
@@ -451,8 +471,8 @@ int	n;
 	double	colv[3], rowv[3];
 	register int	i, j, k;
 
-	if (n < 3 | n > NMBMOD) {
-		fprintf(stderr, "%s: bad number of colors to match\n", progname);
+	if (n < 3) {
+		fprintf(stderr, "%s: too few colors to match!\n", progname);
 		exit(1);
 	}
 	if (n == 3)
@@ -504,18 +524,30 @@ int	n;
 cvtcolor(cout, cin)		/* convert color according to our mapping */
 COLOR	cout, cin;
 {
+	COLOR	ctmp;
+
+	bresp(ctmp, cin);
+	cresp(cout, ctmp);
+	if (colval(cout,RED) < 0.)
+		colval(cout,RED) = 0.;
+	if (colval(cout,GRN) < 0.)
+		colval(cout,GRN) = 0.;
+	if (colval(cout,BLU) < 0.)
+		colval(cout,BLU) = 0.;
+}
+
+
+cresp(cout, cin)		/* transform color according to matrix */
+COLOR	cout, cin;
+{
 	double	r, g, b;
 
-	bresp(cout, cin);
-	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[1][0] + colval(cout,1)*solmat[1][1]
-			+ colval(cout,2)*solmat[1][2];
-	if (g < 0) g = 0;
-	b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1]
-			+ colval(cout,2)*solmat[2][2];
-	if (b < 0) b = 0;
+	r = colval(cin,0)*solmat[0][0] + colval(cin,1)*solmat[0][1]
+			+ colval(cin,2)*solmat[0][2];
+	g = colval(cin,0)*solmat[1][0] + colval(cin,1)*solmat[1][1]
+			+ colval(cin,2)*solmat[1][2];
+	b = colval(cin,0)*solmat[2][0] + colval(cin,1)*solmat[2][1]
+			+ colval(cin,2)*solmat[2][2];
 	setcolor(cout, r, g, b);
 }
 
@@ -537,9 +569,10 @@ register float	xyYin[3];
 
 picdebug()			/* put out debugging picture */
 {
+	static COLOR	blkcol = BLKCOLOR;
 	COLOR	*scan;
-	int	y;
-	register int	x, i;
+	int	y, i;
+	register int	x, rg;
 
 	if (fseek(stdin, 0L, 0) == EOF) {
 		fprintf(stderr, "%s: cannot seek on input picture\n", progname);
@@ -565,11 +598,16 @@ picdebug()			/* put out debugging picture */
 			exit(1);
 		}
 		for (x = 0; x < xmax; x++) {
-			i = chartndx(x, y, CENTCVG);
-			if (i < 0)
+			rg = chartndx(x, y, &i);
+			if (rg == RG_CENT) {
+				if (!(1L<<i & gmtflags) || (x+y)&07)
+					copycolor(scan[x], mbRGB[i]);
+				else
+					copycolor(scan[x], blkcol);
+			} else if (rg == RG_CORR)
 				cvtcolor(scan[x], scan[x]);
-			else
-				copycolor(scan[x], mbRGB[i]);
+			else if (rg != RG_ORIG)
+				copycolor(scan[x], blkcol);
 		}
 		if (fwritescan(scan, xmax, debugfp) < 0) {
 			fprintf(stderr, "%s: error writing debugging picture\n",
@@ -586,16 +624,19 @@ picdebug()			/* put out debugging picture */
 clrdebug()			/* put out debug picture from color input */
 {
 	static COLR	blkclr = BLKCOLR;
-	COLR	mbclr[24], cvclr[24];
+	COLR	mbclr[24], cvclr[24], orclr[24];
 	COLR	*scan;
 	COLOR	ctmp;
-	int	y;
-	register int	i, x;
+	int	y, i;
+	register int	x, rg;
 						/* 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<<i) {
+			setcolr(orclr[i], colval(inpRGB[i],RED),
+					colval(inpRGB[i],GRN),
+					colval(inpRGB[i],BLU));
 			cvtcolor(ctmp, inpRGB[i]);
 			setcolr(cvclr[i], colval(ctmp,RED),
 					colval(ctmp,GRN), colval(ctmp,BLU));
@@ -613,14 +654,20 @@ clrdebug()			/* put out debug picture from color input
 	fprtresolu(xmax, ymax, debugfp);
 						/* write debug picture */
 	for (y = ymax-1; y >= 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<<i)
-				copycolr(scan[x], cvclr[i]);
-			else
+		for (x = 0; x < xmax; x++) {
+			rg = chartndx(x, y, &i);
+			if (rg == RG_CENT) {
+				if (!(1L<<i & gmtflags) || (x+y)&07)
+					copycolr(scan[x], mbclr[i]);
+				else
+					copycolr(scan[x], blkclr);
+			} else if (rg == RG_BORD || !(1L<<i & inpflags))
 				copycolr(scan[x], blkclr);
+			else if (rg == RG_ORIG)
+				copycolr(scan[x], orclr[i]);
+			else /* rg == RG_CORR */
+				copycolr(scan[x], cvclr[i]);
+		}
 		if (fwritecolrs(scan, xmax, debugfp) < 0) {
 			fprintf(stderr, "%s: error writing debugging picture\n",
 					progname);