--- ray/src/common/spec_rgb.c	1996/10/02 11:05:18	2.6
+++ ray/src/common/spec_rgb.c	2003/06/30 19:04:29	2.13
@@ -1,18 +1,33 @@
-/* Copyright (c) 1990 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: spec_rgb.c,v 2.13 2003/06/30 19:04:29 greg Exp $";
 #endif
-
 /*
  * Convert colors and spectral ranges.
+ * Added von Kries white-balance calculations 10/01 (GW).
+ *
+ * Externals declared in color.h
  */
 
+#include "copyright.h"
+
+#include <stdio.h>
+#include <string.h>
 #include "color.h"
 
+#define CEPS	1e-4			/* color epsilon */
 
-RGBPRIMS  stdprims = STDPRIMS;		/* standard primary chromaticities */
+#define CEQ(v1,v2)	((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)
 
+#define XYEQ(c1,c2)	(CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))
+
+
+RGBPRIMS  stdprims = STDPRIMS;	/* standard primary chromaticities */
+
+COLOR  cblack = BLKCOLOR;		/* global black color */
+COLOR  cwhite = WHTCOLOR;		/* global white color */
+
+float  xyneu[2] = {1./3., 1./3.};	/* neutral xy chromaticities */
+
 /*
  *	The following table contains the CIE tristimulus integrals
  *  for X, Y, and Z.  The table is cumulative, so that
@@ -42,7 +57,7 @@ static BYTE  chroma[3][NINC] = {
 	}
 };
 
-COLORMAT  xyz2rgbmat = {		/* XYZ to RGB */
+COLORMAT  xyz2rgbmat = {		/* XYZ to RGB (no white balance) */
 	{(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
 	 (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
 	 (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
@@ -54,7 +69,7 @@ COLORMAT  xyz2rgbmat = {		/* XYZ to RGB */
 	 (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
 };
 
-COLORMAT  rgb2xyzmat = {		/* RGB to XYZ */
+COLORMAT  rgb2xyzmat = {		/* RGB to XYZ (no white balance) */
 	{CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
 	{CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
 	{(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
@@ -62,8 +77,20 @@ COLORMAT  rgb2xyzmat = {		/* RGB to XYZ */
 	 (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
 };
 
+COLORMAT  vkmat = {		/* Sharp primary matrix */
+	{ 1.2694, -0.0988, -0.1706},
+	{-0.8364,  1.8006,  0.0357},
+	{ 0.0297, -0.0315,  1.0018}
+};
 
+COLORMAT  ivkmat = {		/* inverse Sharp primary matrix */
+	{ 0.8156,  0.0472,  0.1372},
+	{ 0.3791,  0.5769,  0.0440},
+	{-0.0123,  0.0167,  0.9955}
+};
 
+
+void
 spec_rgb(col, s, e)		/* compute RGB color from spectral range */
 COLOR  col;
 int  s, e;
@@ -75,6 +102,7 @@ int  s, e;
 }
 
 
+void
 spec_cie(col, s, e)		/* compute a color from a spectral range */
 COLOR  col;		/* returned color */
 int  s, e;		/* starting and ending wavelengths */
@@ -109,23 +137,82 @@ int  s, e;		/* starting and ending wavelengths */
 }
 
 
+void
+cie_rgb(rgb, xyz)		/* convert CIE color to standard RGB */
+COLOR	rgb;
+COLOR  xyz;
+{
+	colortrans(rgb, xyz2rgbmat, xyz);
+	clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
+}
+
+
+int
+clipgamut(col, brt, gamut, lower, upper)	/* clip to gamut cube */
+COLOR  col;
+double  brt;
+int  gamut;
+COLOR  lower, upper;
+{
+	int  rflags = 0;
+	double  brtmin, brtmax, v, vv;
+	COLOR  cgry;
+	register int  i;
+					/* check for no check */
+	if (gamut == 0) return(0);
+					/* check brightness limits */
+	brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
+	if (gamut & CGAMUT_LOWER && brt < brtmin) {
+		copycolor(col, lower);
+		return(CGAMUT_LOWER);
+	}
+	brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
+	if (gamut & CGAMUT_UPPER && brt > brtmax) {
+		copycolor(col, upper);
+		return(CGAMUT_UPPER);
+	}
+					/* compute equivalent grey */
+	v = (brt - brtmin)/(brtmax - brtmin);
+	for (i = 0; i < 3; i++)
+		cgry[i] = v*upper[i] + (1.-v)*lower[i];
+	vv = 1.;			/* check each limit */
+	for (i = 0; i < 3; i++)
+		if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
+			v = (lower[i]+CEPS - cgry[i])/(col[i] - cgry[i]);
+			if (v < vv) vv = v;
+			rflags |= CGAMUT_LOWER;
+		} else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
+			v = (upper[i]-CEPS - cgry[i])/(col[i] - cgry[i]);
+			if (v < vv) vv = v;
+			rflags |= CGAMUT_UPPER;
+		}
+	if (rflags)			/* desaturate to cube face */
+		for (i = 0; i < 3; i++)
+			col[i] = vv*col[i] + (1.-vv)*cgry[i];
+	return(rflags);
+}
+
+
+void
 colortrans(c2, mat, c1)		/* convert c1 by mat and put into c2 */
+register COLOR  c2;
 register COLORMAT  mat;
-register COLOR  c1, c2;
+register COLOR  c1;
 {
-	static float  cout[3];
+	COLOR	cout;
 
 	cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
 	cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
 	cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
-	if((c2[0] = cout[0]) < 0.) c2[0] = 0.;
-	if((c2[1] = cout[1]) < 0.) c2[1] = 0.;
-	if((c2[2] = cout[2]) < 0.) c2[2] = 0.;
+
+	copycolor(c2, cout);
 }
 
 
+void
 multcolormat(m3, m2, m1)	/* multiply m1 by m2 and put into m3 */
-COLORMAT  m1, m2, m3;		/* m3 can be either m1 or m2 w/o harm */
+COLORMAT  m3;			/* m3 can be either m1 or m2 w/o harm */
+COLORMAT  m2, m1;
 {
 	COLORMAT  mt;
 	register int  i, j;
@@ -139,6 +226,7 @@ COLORMAT  m1, m2, m3;		/* m3 can be either m1 or m2 w/
 }
 
 
+void
 compxyz2rgbmat(mat, pr)		/* compute conversion from CIE to RGB space */
 COLORMAT  mat;
 register RGBPRIMS  pr;
@@ -189,6 +277,7 @@ register RGBPRIMS  pr;
 }
 
 
+void
 comprgb2xyzmat(mat, pr)		/* compute conversion from RGB to CIE space */
 COLORMAT  mat;
 register RGBPRIMS  pr;
@@ -226,6 +315,7 @@ register RGBPRIMS  pr;
 }
 
 
+void
 comprgb2rgbmat(mat, pr1, pr2)	/* compute conversion from RGB1 to RGB2 */
 COLORMAT  mat;
 RGBPRIMS  pr1, pr2;
@@ -242,4 +332,84 @@ RGBPRIMS  pr1, pr2;
 	compxyz2rgbmat(xyztopr2, pr2);
 				/* combine transforms */
 	multcolormat(mat, pr1toxyz, xyztopr2);
+}
+
+
+void
+compxyzWBmat(mat, wht1, wht2)	/* CIE von Kries transform from wht1 to wht2 */
+COLORMAT  mat;
+float  wht1[2], wht2[2];
+{
+	COLOR	cw1, cw2;
+	if (XYEQ(wht1,wht2)) {
+		mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
+		mat[0][1] = mat[0][2] = mat[1][0] =
+		mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
+		return;
+	}
+	cw1[RED] = wht1[CIEX]/wht1[CIEY];
+	cw1[GRN] = 1.;
+	cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
+	colortrans(cw1, vkmat, cw1);
+	cw2[RED] = wht2[CIEX]/wht2[CIEY];
+	cw2[GRN] = 1.;
+	cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
+	colortrans(cw2, vkmat, cw2);
+	mat[0][0] = cw2[RED]/cw1[RED];
+	mat[1][1] = cw2[GRN]/cw1[GRN];
+	mat[2][2] = cw2[BLU]/cw1[BLU];
+	mat[0][1] = mat[0][2] = mat[1][0] =
+	mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
+	multcolormat(mat, vkmat, mat);
+	multcolormat(mat, mat, ivkmat);
+}
+
+
+void
+compxyz2rgbWBmat(mat, pr)	/* von Kries conversion from CIE to RGB space */
+COLORMAT  mat;
+RGBPRIMS  pr;
+{
+	COLORMAT	wbmat;
+
+	compxyz2rgbmat(mat, pr);
+	if (XYEQ(pr[WHT],xyneu))
+		return;
+	compxyzWBmat(wbmat, xyneu, pr[WHT]);
+	multcolormat(mat, wbmat, mat);
+}
+
+void
+comprgb2xyzWBmat(mat, pr)	/* von Kries conversion from RGB to CIE space */
+COLORMAT  mat;
+RGBPRIMS  pr;
+{
+	COLORMAT	wbmat;
+	
+	comprgb2xyzmat(mat, pr);
+	if (XYEQ(pr[WHT],xyneu))
+		return;
+	compxyzWBmat(wbmat, pr[WHT], xyneu);
+	multcolormat(mat, mat, wbmat);
+}
+
+void
+comprgb2rgbWBmat(mat, pr1, pr2)	/* von Kries conversion from RGB1 to RGB2 */
+COLORMAT  mat;
+RGBPRIMS  pr1, pr2;
+{
+	COLORMAT  pr1toxyz, xyztopr2, wbmat;
+
+	if (pr1 == pr2) {
+		mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
+		mat[0][1] = mat[0][2] = mat[1][0] =
+		mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
+		return;
+	}
+	comprgb2xyzmat(pr1toxyz, pr1);
+	compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]);
+	compxyz2rgbmat(xyztopr2, pr2);
+				/* combine transforms */
+	multcolormat(mat, pr1toxyz, wbmat);
+	multcolormat(mat, mat, xyztopr2);
 }