--- ray/src/px/clrtab.c	1992/10/12 12:58:39	2.1
+++ ray/src/px/clrtab.c	2005/09/19 02:23:58	2.18
@@ -1,16 +1,18 @@
-/* Copyright (c) 1992 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: clrtab.c,v 2.18 2005/09/19 02:23:58 greg Exp $";
 #endif
-
 /*
  * Simple median-cut color quantization based on colortab.c
  */
 
-#include "standard.h"
+#include "copyright.h"
 
+#include <string.h>
+
+#include "standard.h"
 #include "color.h"
+#include "clrtab.h"
+
 				/* histogram resolution */
 #define NRED		36
 #define NGRN		48
@@ -23,32 +25,55 @@ static char SCCSid[] = "$SunId$ LBL";
 #define part(cn)	((cn)>>2)
 #define prim(cn)	((cn)&3)
 				/* our color table (global) */
-BYTE	clrtab[256][3];
+extern BYTE	clrtab[256][3];
 				/* histogram of colors / color assignments */
 static unsigned	histo[NRED][NGRN][NBLU];
 #define cndx(c)		histo[((c)[RED]*NRED)>>8][((c)[GRN]*NGRN)>>8][((c)[BLU]*NBLU)>>8]
 				/* initial color cube boundary */
-static int	CLRCUBE[3][2] = {0,NRED,0,NGRN,0,NBLU};
+static int	CLRCUBE[3][2] = {{0,NRED},{0,NGRN},{0,NBLU}};
 				/* maximum propagated error during dithering */
 #define MAXERR		20
+				/* define CLOSEST to get closest colors */
+#ifndef CLOSEST
+#define CLOSEST		1	/* this step takes a little longer */
+#endif
 
 
-new_histo()		/* clear our histogram */
+#ifdef CLOSEST
+static void closest(int n);
+static void setclosest(BYTE *nl[], int r, int g, int b);
+static void addneigh(BYTE *nl[], int i, int j);
+static unsigned int dist(BYTE col[3], int r, int g, int b);
+#endif
+static void cut(int box[3][2], int c0, int c1);
+static int split(int box[3][2]);
+static void mktabent(int p, int box[3][2]);
+
+
+extern int
+new_histo(		/* clear our histogram */
+	int	n
+)
 {
-	bzero((char *)histo, sizeof(histo));
+	memset((void *)histo, '\0', sizeof(histo));
+	return(0);
 }
 
 
-cnt_pixel(col)		/* add pixel to our histogram */
-register BYTE	col[];
+extern void
+cnt_pixel(		/* add pixel to our histogram */
+	register BYTE	col[]
+)
 {
 	cndx(col)++;
 }
 
 
-cnt_colrs(cs, n)		/* add a scanline to our histogram */
-register COLR	*cs;
-register int	n;
+extern void
+cnt_colrs(		/* add a scanline to our histogram */
+	register COLR	*cs,
+	register int	n
+)
 {
 	while (n-- > 0) {
 		cndx(cs[0])++;
@@ -57,8 +82,10 @@ register int	n;
 }
 
 
-new_clrtab(ncolors)		/* make new color table using ncolors */
-int	ncolors;
+extern int
+new_clrtab(		/* make new color table using ncolors */
+	int	ncolors
+)
 {
 	if (ncolors < 1)
 		return(0);
@@ -66,23 +93,31 @@ int	ncolors;
 		ncolors = 256;
 				/* partition color space */
 	cut(CLRCUBE, 0, ncolors);
+#ifdef CLOSEST
+	closest(ncolors);	/* ensure colors picked are closest */
+#endif
+				/* reset dithering function */
+	dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
 				/* return new color table size */
 	return(ncolors);
 }
 
 
-int
-map_pixel(col)			/* get pixel for color */
-register BYTE	col[];
+extern int
+map_pixel(			/* get pixel for color */
+	register BYTE	col[]
+)
 {
     return(cndx(col));
 }
 
 
-map_colrs(bs, cs, n)		/* convert a scanline to color index values */
-register BYTE	*bs;
-register COLR	*cs;
-register int	n;
+extern void
+map_colrs(		/* convert a scanline to color index values */
+	register BYTE	*bs,
+	register COLR	*cs,
+	register int	n
+)
 {
 	while (n-- > 0) {
 		*bs++ = cndx(cs[0]);
@@ -91,21 +126,24 @@ register int	n;
 }
 
 
-dith_colrs(bs, cs, n)		/* convert scanline to dithered index values */
-register BYTE	*bs;
-register COLR	*cs;
-int	n;
+extern void
+dith_colrs(		/* convert scanline to dithered index values */
+	register BYTE	*bs,
+	register COLR	*cs,
+	int	n
+)
 {
-	static short	(*cerr)[3];
+	static short	(*cerr)[3] = NULL;
 	static int	N = 0;
 	int	err[3], errp[3];
 	register int	x, i;
 
 	if (n != N) {		/* get error propogation array */
-		if (N)
-			cerr = (short (*)[3])realloc((char *)cerr,
-					3*n*sizeof(short));
-		else
+		if (N) {
+			free((void *)cerr);
+			cerr = NULL;
+		}
+		if (n)
 			cerr = (short (*)[3])malloc(3*n*sizeof(short));
 		if (cerr == NULL) {
 			N = 0;
@@ -113,7 +151,7 @@ int	n;
 			return;
 		}
 		N = n;
-		bzero((char *)cerr, 3*N*sizeof(short));
+		memset((void *)cerr, '\0', 3*N*sizeof(short));
 	}
 	err[0] = err[1] = err[2] = 0;
 	for (x = 0; x < n; x++) {
@@ -138,10 +176,12 @@ int	n;
 }
 
 
-static
-cut(box, c0, c1)			/* partition color space */
-register int	box[3][2];
-int	c0, c1;
+static void
+cut(			/* partition color space */
+	register int	box[3][2],
+	int	c0,
+	int	c1
+)
 {
 	register int	branch;
 	int	kb[3][2];
@@ -152,7 +192,7 @@ int	c0, c1;
 	}
 					/* split box */
 	branch = split(box);
-	bcopy((char *)box, (char *)kb, sizeof(kb));
+	memcpy((void *)kb, (void *)box, sizeof(kb));
 						/* do left (lesser) branch */
 	kb[prim(branch)][1] = part(branch);
 	cut(kb, c0, (c0+c1)>>1);
@@ -164,13 +204,14 @@ int	c0, c1;
 
 
 static int
-split(box)				/* find median cut for box */
-register int	box[3][2];
+split(				/* find median cut for box */
+	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])
@@ -221,21 +262,23 @@ register int	box[3][2];
 }
 
 
-static
-mktabent(p, box)	/* compute average color for box and assign */
-int	p;
-register int	box[3][2];
+static void
+mktabent(	/* compute average color for box and assign */
+	int	p,
+	register int	box[3][2]
+)
 {
-	long	sum[3];
-	int	r, g, n;
-	register int	b, c;
+	unsigned long	sum[3];
+	unsigned	r, g;
+	unsigned long	n;
+	register unsigned	b, c;
 						/* sum pixels in box */
 	n = 0;
 	sum[RED] = sum[GRN] = sum[BLU] = 0;
 	for (r = box[RED][0]; r < box[RED][1]; r++)
 	    for (g = box[GRN][0]; g < box[GRN][1]; g++)
 	    	for (b = box[BLU][0]; b < box[BLU][1]; b++) {
-	    	    if (c = histo[r][g][b]) {
+	    	    if ( (c = histo[r][g][b]) ) {
 			n += c;
 			sum[RED] += (long)c*r;
 			sum[GRN] += (long)c*g;
@@ -243,6 +286,12 @@ register int	box[3][2];
 		    }
 		    histo[r][g][b] = p;		/* assign pixel */
 		}
+	if (n >= (1L<<23)/HMAX) {		/* avoid overflow */
+		sum[RED] /= n;
+		sum[GRN] /= n;
+		sum[BLU] /= n;
+		n = 1;
+	}
 	if (n) {				/* compute average */
 		clrtab[p][RED] = sum[RED]*256/NRED/n;
 		clrtab[p][GRN] = sum[GRN]*256/NGRN/n;
@@ -253,3 +302,121 @@ register int	box[3][2];
 		clrtab[p][BLU] = (box[BLU][0]+box[BLU][1])*256/NBLU/2;
 	}
 }
+
+
+#ifdef CLOSEST
+#define NBSIZ		32
+static void
+closest(			/* make sure we have the closest colors */
+	int	n
+)
+{
+	BYTE	*neigh[256];
+	register int	r, g, b;
+#define i r
+					/* get space for neighbor lists */
+	for (i = 0; i < n; i++) {
+		if ((neigh[i] = (BYTE *)malloc(NBSIZ)) == NULL) {
+			while (i--)
+				free(neigh[i]);
+			return;			/* ENOMEM -- abandon effort */
+		}
+		neigh[i][0] = i;		/* identity is terminator */
+	}
+					/* make neighbor lists */
+	for (r = 0; r < NRED; r++)
+	    for (g = 0; g < NGRN; g++)
+		for (b = 0; b < NBLU; b++) {
+		    if (r < NRED-1 && histo[r][g][b] != histo[r+1][g][b])
+			addneigh(neigh, histo[r][g][b], histo[r+1][g][b]);
+		    if (g < NGRN-1 && histo[r][g][b] != histo[r][g+1][b])
+			addneigh(neigh, histo[r][g][b], histo[r][g+1][b]);
+		    if (b < NBLU-1 && histo[r][g][b] != histo[r][g][b+1])
+			addneigh(neigh, histo[r][g][b], histo[r][g][b+1]);
+		}
+					/* assign closest values */
+	for (r = 0; r < NRED; r++)
+	    for (g = 0; g < NGRN; g++)
+		for (b = 0; b < NBLU; b++)
+		    setclosest(neigh, r, g, b);
+					/* free neighbor lists */
+	for (i = 0; i < n; i++)
+		free(neigh[i]);
+#undef i
+}
+
+
+static void
+addneigh(		/* i and j are neighbors; add them to list */
+	register BYTE	*nl[],
+	register int	i,
+	int	j
+)
+{
+	int	nc;
+	char	*nnl;
+	register int	t;
+	
+	for (nc = 0; nc < 2; nc++) {		/* do both neighbors */
+		for (t = 0; nl[i][t] != i; t++)
+			if (nl[i][t] == j)
+				break;		/* in list already */
+		if (nl[i][t] == i) {		/* add to list */
+			nl[i][t++] = j;
+			if (t % NBSIZ == 0) {	/* enlarge list */
+				if ((nnl = realloc((void *)nl[i],
+						t+NBSIZ)) == NULL)
+					t--;
+				else
+					nl[i] = (BYTE *)nnl;
+			}
+			nl[i][t] = i;		/* terminator */
+		}
+		t = i; i = j; j = t;		/* swap and do it again */
+	}
+}
+
+
+static unsigned int
+dist(		/* find distance from clrtab entry to r,g,b */
+	register BYTE	col[3],
+	int	r,
+	int	g,
+	int	b
+)
+{
+	register int	tmp;
+	register unsigned	sum;
+	
+	tmp = col[RED]*NRED/256 - r;
+	sum = tmp*tmp;
+	tmp = col[GRN]*NGRN/256 - g;
+	sum += tmp*tmp;
+	tmp = col[BLU]*NBLU/256 - b;
+	sum += tmp*tmp;
+	return(sum);
+}
+
+
+static void
+setclosest(		/* find index closest to color and assign */
+	BYTE	*nl[],
+	int	r,
+	int	g,
+	int	b
+)
+{
+	int	ident;
+	unsigned	min;
+	register unsigned	d;
+	register BYTE	*p;
+					/* get starting value */
+	min = dist(clrtab[ident=histo[r][g][b]], r, g, b);
+					/* find minimum */
+	for (p = nl[ident]; *p != ident; p++)
+		if ((d = dist(clrtab[*p], r, g, b)) < min) {
+			min = d;
+			histo[r][g][b] = *p;
+		}
+}
+#endif