--- ray/src/px/macbethcal.c	1995/10/11 10:40:13	2.1
+++ ray/src/px/macbethcal.c	1995/10/17 11:05:17	2.6
@@ -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,31 @@ 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,3,2};
-#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		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 NMBSAT		8	/* Number of MacBeth saturated colors */
+short	mbsat[NMBSAT] = {
+		Red,Green,Blue,Magenta,Yellow,Cyan,
+		Orange,Purple
+	};
+
 int	xmax, ymax;		/* input image dimensions */
 int	bounds[4][2];		/* image coordinates of chart corners */
 double	imgxfm[3][3];		/* coordinate transformation matrix */
 
 COLOR	picRGB[24];		/* picture colors */
 
-double	bramp[NMBNEU][2];	/* brightness ramp */
+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();
@@ -125,6 +159,10 @@ 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);
@@ -138,9 +176,6 @@ userr:
 init()				/* initialize */
 {
 	double	quad[4][2];
-	COLOR	ctmp;
-	double	d;
-	int	i;
 					/* make coordinate transformation */
 	quad[0][0] = bounds[0][0];
 	quad[0][1] = bounds[0][1];
@@ -155,14 +190,6 @@ init()				/* initialize */
 		fprintf(stderr, "%s: bad chart boundaries\n", progname);
 		exit(1);
 	}
-					/* 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);
-	}
 }
 
 
@@ -238,79 +265,123 @@ getcolors()				/* load in picture colors */
 }
 
 
-double
-bresp(x)		/* piecewise linear interpolation of brightness */
-double	x;
+bresp(y, x)		/* piecewise linear interpolation of primaries */
+COLOR	y, x;
 {
-	register int	n = NMBNEU;
+	double	cv[3];
+	register int	i, n;
 
-	while (n > 0 && x < bramp[--n][0])
-		;
-	return( ((bramp[n+1][0] - x)*bramp[n][1] +
-				(x - bramp[n][0])*bramp[n+1][1]) /
-			(bramp[n+1][0] - bramp[n][0]) );
+	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(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]);
 }
 
 
-putmapping()			/* compute and print mapping for pcomb -f */
+compute()			/* compute color mapping */
 {
-	float	clrin[NMBMOD][3], clrout[NMBMOD][3];
-	register int	i, j;
+	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);
+	}
 					/* compute piecewise luminance curve */
 	for (i = 0; i < NMBNEU; i++) {
-		bramp[i][0] = bright(picRGB[mbneu[i]]);
-		bramp[i][1] = bright(mbRGB[mbneu[i]]);
+		copycolor(bramp[i][0], picRGB[mbneu[i]]);
+		copycolor(bramp[i][1], mbRGB[mbneu[i]]);
 	}
 					/* compute color matrix */
-	for (i = 0; i < NMBMOD; i++)
-		for (j = 0; j < 3; j++) {
-			clrin[i][j] = bresp(picRGB[mbmod[i]][j]);
-			clrout[i][j] = mbRGB[mbmod[i]][j];
-		}
+	for (i = 0; i < NMBMOD; i++) {
+		bresp(clrin[i], picRGB[mbmod[i]]);
+		copycolor(clrout[i], mbRGB[mbmod[i]]);
+	}
 	compsoln(clrin, clrout, NMBMOD);
+}
+
+
+putmapping()			/* put out color mapping for pcomb -f */
+{
+	static char	cchar[3] = {'r', 'g', 'b'};
+	register int	i, j;
 					/* print brightness mapping */
-	printf("xval(i) : select(i, 0");
-	for (i = 0; i < NMBNEU; i++)
-		printf(", %g", bramp[i][0]);
-	printf(");\n");
-	printf("yval(i) : select(i, 0");
-	for (i = 0; i < NMBNEU; i++)
-		printf(", %g", bramp[i][1]);
-	printf(");\n");
-	printf("ifind(x,f,n) : if(1.5-n, 1, if(x-f(n), n, ifind(x,f,n-1)));\n");
-	printf("binterp(i,x) : ((xval(i+1)-x)*yval(i) + (x-xval(i))*yval(i+1))/\n");
-	printf("\t\t(xval(i+1) - xval(i));\n");
-	printf("bresp(x) : binterp(ifind(x,xval,xval(0)-1), x);\n");
-	printf("nred = bresp(ri(1));\n");
-	printf("ngrn = bresp(gi(1));\n");
-	printf("nblu = bresp(bi(1));\n");
+	for (j = 0; j < 3; j++) {
+		printf("%cxa(i) : select(i", cchar[j]);
+		for (i = 0; i < NMBNEU; i++)
+			printf(",%g", colval(bramp[i][0],j));
+		printf(");\n");
+		printf("%cya(i) : select(i", cchar[j]);
+		for (i = 0; i < NMBNEU; i++)
+			printf(",%g", colval(bramp[i][1],j));
+		printf(");\n");
+		printf("%c = %ci(1);\n", cchar[j], cchar[j]);
+		printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
+				cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
+				cchar[j], cchar[j]);
+		printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
+		printf("%cn = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
+				cchar[j], cchar[j], cchar[j],
+				cchar[j], cchar[j], cchar[j]);
+		printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
+				cchar[j], cchar[j], cchar[j],
+				cchar[j], cchar[j]);
+		printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
+				cchar[j], cchar[j], cchar[j], cchar[j]);
+	}
 					/* print color mapping */
-	printf("ro = %g*nred + %g*ngrn + %g*nblu\n",
-			solmat[0][0], solmat[1][0], solmat[2][0]);
-	printf("go = %g*nred + %g*ngrn + %g*nblu\n",
-			solmat[0][1], solmat[1][1], solmat[2][1]);
-	printf("bo = %g*nred + %g*ngrn + %g*nblu\n",
-			solmat[0][2], solmat[1][2], solmat[2][2]);
+	printf("ro = %g*rn + %g*gn + %g*bn ;\n",
+			solmat[0][0], solmat[0][1], solmat[0][2]);
+	printf("go = %g*rn + %g*gn + %g*bn ;\n",
+			solmat[1][0], solmat[1][1], solmat[1][2]);
+	printf("bo = %g*rn + %g*gn + %g*bn ;\n",
+			solmat[2][0], solmat[2][1], solmat[2][2]);
 }
 
 
-compsoln(cin, cout, n)		/* solve 3x3 system */
-float	cin[][3], cout[][3];
+compsoln(cin, cout, n)		/* solve 3xN system using least-squares */
+COLOR	cin[], cout[];
 int	n;
 {
 	extern double	mx3d_adjoint(), fabs();
 	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] = 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",
@@ -324,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] = cout[j][i];
-		mx3d_transform(rowv, invmat, colv);
-		for (j = 0; j < 3; j++)
-			solmat[j][i] = colv[j];
+			solmat[i][j] = rowv[j];
 	}
 }
 
@@ -337,18 +416,18 @@ cvtcolor(cout, cin)		/* convert color according to our
 COLOR	cout, cin;
 {
 	double	r, g, b;
-	double	r1, g1, b1;
 
-	r = bresp(colval(cin,RED));
-	g = bresp(colval(cin,GRN));
-	b = bresp(colval(cin,BLU));
-	r1 = r*solmat[0][0] + g*solmat[1][0] + b*solmat[2][0];
-	if (r1 < 0) r1 = 0;
-	g1 = r*solmat[0][1] + g*solmat[1][1] + b*solmat[2][1];
-	if (g1 < 0) g1 = 0;
-	b1 = r*solmat[0][2] + g*solmat[1][2] + b*solmat[2][2];
-	if (b1 < 0) b1 = 0;
-	setcolor(cout, r1, g1, b1);
+	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;
+	setcolor(cout, r, g, b);
 }
 
 
@@ -373,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--) {