--- ray/src/rt/colortab.c	1989/10/02 17:12:42	1.2
+++ ray/src/rt/colortab.c	2003/02/22 02:07:28	2.5
@@ -1,9 +1,6 @@
-/* Copyright (c) 1989 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: colortab.c,v 2.5 2003/02/22 02:07:28 greg Exp $";
 #endif
-
 /*
  * colortab.c - allocate and control dynamic color table.
  *
@@ -14,140 +11,251 @@ static char SCCSid[] = "$SunId$ LBL";
  *	histogram is cleared.  This algorithm
  *	performs only as well as the next drawing's color
  *	distribution is correlated to the last.
+ *
+ *  External symbols declared in drvier.h
  */
 
-#include "color.h"
+/* ====================================================================
+ * The Radiance Software License, Version 1.0
+ *
+ * Copyright (c) 1990 - 2002 The Regents of the University of California,
+ * through Lawrence Berkeley National Laboratory.   All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ *         notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in
+ *       the documentation and/or other materials provided with the
+ *       distribution.
+ *
+ * 3. The end-user documentation included with the redistribution,
+ *           if any, must include the following acknowledgment:
+ *             "This product includes Radiance software
+ *                 (http://radsite.lbl.gov/)
+ *                 developed by the Lawrence Berkeley National Laboratory
+ *               (http://www.lbl.gov/)."
+ *       Alternately, this acknowledgment may appear in the software itself,
+ *       if and wherever such third-party acknowledgments normally appear.
+ *
+ * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
+ *       and "The Regents of the University of California" must
+ *       not be used to endorse or promote products derived from this
+ *       software without prior written permission. For written
+ *       permission, please contact radiance@radsite.lbl.gov.
+ *
+ * 5. Products derived from this software may not be called "Radiance",
+ *       nor may "Radiance" appear in their name, without prior written
+ *       permission of Lawrence Berkeley National Laboratory.
+ *
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
+ * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+ * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ * ====================================================================
+ *
+ * This software consists of voluntary contributions made by many
+ * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
+ * information on Lawrence Berkeley National Laboratory, please see
+ * <http://www.lbl.gov/>.
+ */
 
-#define NULL		0
+#include "standard.h"
+
+#include "color.h"
 				/* histogram resolution */
-#define NRED		28
+#define NRED		24
 #define NGRN		32
-#define NBLU		22
+#define NBLU		16
 #define HMAX		NGRN
 				/* minimum box count for adaptive partition */
 #define MINSAMP		7
-				/* maximum frame buffer depth */
-#define FBDEPTH		8
-				/* color map resolution */
-#define MAPSIZ		128
+				/* maximum distance^2 before color reassign */
+#define MAXDST2		12
 				/* map a color */
 #define map_col(c,p)	clrmap[p][ colval(c,p)<1. ? \
-				(int)(colval(c,p)*MAPSIZ) : MAPSIZ-1 ]
+				(int)(colval(c,p)*256.) : 255 ]
 				/* color partition tree */
 #define CNODE		short
-#define set_branch(p,c)	((c)<<2|(p))
+#define set_branch(p,c) ((c)<<2|(p))
 #define set_pval(pv)	((pv)<<2|3)
+#define is_branch(cn)	(((cn)&3)!=3)
 #define is_pval(cn)	(((cn)&3)==3)
 #define part(cn)	((cn)>>2)
 #define prim(cn)	((cn)&3)
 #define pval(cn)	((cn)>>2)
 				/* our color table */
-COLR	clrtab[1<<FBDEPTH];
+static struct tabent {
+	long	sum[3];		/* sum of colors using this entry */
+	int	n;		/* number of colors */
+	BYTE	ent[3];		/* current table value */
+}	*clrtab = NULL;
+				/* color cube partition */
+static CNODE	*ctree = NULL;
 				/* our color correction map */
-static BYTE	clrmap[3][MAPSIZ];
+static BYTE	clrmap[3][256];
 				/* histogram of colors used */
-static unsigned	histo[NRED][NGRN][NBLU];
-				/* initial color cube boundaries */
-static int	CLRCUBE[3][2] = {0,NRED,0,NGRN,0,NBLU};
-				/* color cube partition */
-static CNODE	ctree[1<<(FBDEPTH+1)];
+static unsigned short	histo[NRED][NGRN][NBLU];
+				/* initial color cube boundary */
+static int	CLRCUBE[3][2] = {{0,NRED},{0,NGRN},{0,NBLU}};
 
+static int	split();
+static void	cut();
 
-COLR *
-get_ctab(ncolors)	/* assign a color table with max ncolors */
+
+int
+new_ctab(ncolors)		/* start new color table with max ncolors */
 int	ncolors;
 {
-	if (ncolors < 1 || ncolors > 1<<FBDEPTH)
-		return(NULL);
-				/* partition color table */
-	cut(ctree, FBDEPTH, CLRCUBE, 0, ncolors);
+	int	treesize;
+
+	if (ncolors < 1)
+		return(0);
+				/* free old tables */
+	if (clrtab != NULL)
+		free((void *)clrtab);
+	if (ctree != NULL)
+		free((void *)ctree);
+				/* get new tables */
+	for (treesize = 1; treesize < ncolors; treesize <<= 1)
+		;
+	treesize <<= 1;
+	clrtab = (struct tabent *)calloc(ncolors, sizeof(struct tabent));
+	ctree = (CNODE *)malloc(treesize*sizeof(CNODE));
+	if (clrtab == NULL || ctree == NULL)
+		return(0);
+				/* partition color space */
+	cut(ctree, 0, CLRCUBE, 0, ncolors);
 				/* clear histogram */
-	bzero(histo, sizeof(histo));
-				/* return color table */
-	return(clrtab);
+	bzero((char *)histo, sizeof(histo));
+				/* return number of colors used */
+	return(ncolors);
 }
 
 
 int
-get_pixel(col)			/* get pixel for color */
+get_pixel(col, set_pixel)	/* get pixel for color */
 COLOR	col;
+void	(*set_pixel)();
 {
+	int	r, g, b;
 	int	cv[3];
 	register CNODE	*tp;
 	register int	h;
-					/* map color */
-	cv[RED] = map_col(col,RED);
-	cv[GRN] = map_col(col,GRN);
-	cv[BLU] = map_col(col,BLU);
-					/* add to histogram */
+						/* map color */
+	r = map_col(col,RED);
+	g = map_col(col,GRN);
+	b = map_col(col,BLU);
+						/* reduce resolution */
+	cv[RED] = (r*NRED)>>8;
+	cv[GRN] = (g*NGRN)>>8;
+	cv[BLU] = (b*NBLU)>>8;
+						/* add to histogram */
 	histo[cv[RED]][cv[GRN]][cv[BLU]]++;
-					/* find pixel value in tree */
-	tp = ctree;
-	for (h = FBDEPTH; h > 0; h--) {
-		if (is_pval(*tp))
-			break;
+						/* find pixel in tree */
+	for (tp = ctree, h = 0; is_branch(*tp); h++)
 		if (cv[prim(*tp)] < part(*tp))
-			tp++;			/* left branch */
+			tp += 1<<h;		/* left branch */
 		else
-			tp += 1<<h;		/* right branch */
-	}
-#ifdef notdef
-	printf("distance (%d,%d,%d)\n",
-		cv[RED] - clrtab[pval(*tp)][RED]*NRED/256, 
-		cv[GRN] - clrtab[pval(*tp)][GRN]*NGRN/256,
-		cv[BLU] - clrtab[pval(*tp)][BLU]*NBLU/256);
+			tp += 1<<(h+1);		/* right branch */
+	h = pval(*tp);
+						/* add to color table */
+	clrtab[h].sum[RED] += r;
+	clrtab[h].sum[GRN] += g;
+	clrtab[h].sum[BLU] += b;
+	clrtab[h].n++;
+					/* recompute average */
+	r = clrtab[h].sum[RED] / clrtab[h].n;
+	g = clrtab[h].sum[GRN] / clrtab[h].n;
+	b = clrtab[h].sum[BLU] / clrtab[h].n;
+					/* check for movement */
+	if (clrtab[h].n == 1 ||
+			(r-clrtab[h].ent[RED])*(r-clrtab[h].ent[RED]) +
+			(g-clrtab[h].ent[GRN])*(g-clrtab[h].ent[GRN]) +
+			(b-clrtab[h].ent[BLU])*(b-clrtab[h].ent[BLU]) > MAXDST2) {
+		clrtab[h].ent[RED] = r;
+		clrtab[h].ent[GRN] = g; /* reassign pixel */
+		clrtab[h].ent[BLU] = b;
+#ifdef DEBUG
+		sprintf(errmsg, "pixel %d = (%d,%d,%d) (%d refs)\n",
+				h, r, g, b, clrtab[h].n);
+		eputs(errmsg);
 #endif
-	return(pval(*tp));
+		(*set_pixel)(h, r, g, b);
+	}
+	return(h);				/* return pixel value */
 }
 
 
-make_cmap(gam)			/* make gamma correction map */
-double  gam;
+void
+make_gmap(gam)			/* make gamma correction map */
+double	gam;
 {
-	extern double	pow();
-	double	d;
 	register int	i;
 	
-	for (i = 0; i < MAPSIZ; i++) {
-		d = pow(i/(double)MAPSIZ, 1.0/gam);
-		clrmap[RED][i] = d * NRED;
-		clrmap[GRN][i] = d * NGRN;
-		clrmap[BLU][i] = d * NBLU;
-	}
+	for (i = 0; i < 256; i++)
+		clrmap[RED][i] =
+		clrmap[GRN][i] =
+		clrmap[BLU][i] = 256.0 * pow((i+0.5)/256.0, 1.0/gam);
 }
 
 
-static
-cut(tree, height, box, c0, c1)		/* partition color space */
+void
+set_cmap(rmap, gmap, bmap)	/* set custom color correction map */
+BYTE	*rmap, *gmap, *bmap;
+{
+	bcopy((char *)rmap, (char *)clrmap[RED], 256);
+	bcopy((char *)gmap, (char *)clrmap[GRN], 256);
+	bcopy((char *)bmap, (char *)clrmap[BLU], 256);
+}
+
+
+void
+map_color(rgb, col)		/* map a color to a byte triplet */
+BYTE	rgb[3];
+COLOR	col;
+{
+	rgb[RED] = map_col(col,RED);
+	rgb[GRN] = map_col(col,GRN);
+	rgb[BLU] = map_col(col,BLU);
+}
+
+
+static void
+cut(tree, level, box, c0, c1)		/* partition color space */
 register CNODE	*tree;
-int	height;
+int	level;
 register int	box[3][2];
 int	c0, c1;
 {
 	int	kb[3][2];
 	
-	if (c1-c0 == 1) {		/* assign color */
-		clrtab[c0][RED] = ((box[RED][0]+box[RED][1])<<7)/NRED;
-		clrtab[c0][GRN] = ((box[GRN][0]+box[GRN][1])<<7)/NGRN;
-		clrtab[c0][BLU] = ((box[BLU][0]+box[BLU][1])<<7)/NBLU;
-		clrtab[c0][EXP] = COLXS;
+	if (c1-c0 <= 1) {		/* assign pixel */
 		*tree = set_pval(c0);
-#ifdef notdef
-		printf("final box size = (%d,%d,%d)\n",
-			box[RED][1] - box[RED][0],
-			box[GRN][1] - box[GRN][0],
-			box[BLU][1] - box[BLU][0]);
-#endif
 		return;
 	}
 					/* split box */
 	*tree = split(box);
-	bcopy(box, kb, sizeof(kb));
+	bcopy((char *)box, (char *)kb, sizeof(kb));
+						/* do left (lesser) branch */
 	kb[prim(*tree)][1] = part(*tree);
-	cut(tree+1, height-1, kb, c0, (c0+c1)>>1);		/* lesser */
+	cut(tree+(1<<level), level+1, kb, c0, (c0+c1)>>1);
+						/* do right branch */
 	kb[prim(*tree)][0] = part(*tree);
 	kb[prim(*tree)][1] = box[prim(*tree)][1];
-	cut(tree+(1<<height), height-1, kb, (c0+c1)>>1, c1);	/* greater */
+	cut(tree+(1<<(level+1)), level+1, kb, (c0+c1)>>1, c1);
 }
 
 
@@ -158,7 +266,7 @@ register int	box[3][2];
 #define c0	r
 	register int	r, g, b;
 	int	pri;
-	int	t[HMAX], med;
+	long	t[HMAX], med;
 					/* find dominant axis */
 	pri = RED;
 	if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])