[comp.graphics] Son of 'Variance-based' Color Quantization

craig@weedeater.math.yale.edu (Craig Kolb) (08/26/89)

As it's starting to look as if the color quantization discussion is in
danger of rearing its head again, I thought I'd take this opportunity
to post my latest implementation of the Wan/Wong/Prusinkiewicz variation
on Heckbert's algorithm.  Rather than a huge, ugly program, this time around
the algorithm is implemented as a routine which, given a prequantized
'full color' image, computes a colormap and an RGB-->index map and returns
the length of the computed colormap.

Included are a hastily-written manual page and a demo program which reads a
Utah Raster format RLE image, quantizes it, and displays the
results on an Iris 4D series workstation.

As always, bug fixes and/or improvements are greatly appreciated.

Craig Kolb
kolb@yale.edu

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of archive 1 (of 1)."
# Contents:  Makefile PORTING colorquant.3 colorquant.c quantdemo.c
# Wrapped by craig@weedeater on Fri Aug 25 15:50:37 1989
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Makefile' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Makefile'\"
else
echo shar: Extracting \"'Makefile'\" \(711 characters\)
sed "s/^X//" >'Makefile' <<'END_OF_FILE'
X#
X# Makefile for variance-based color quantization routine
X#
X# Craig Kolb 12/88
X#
X# The demo program, 'quantdemo', uses the Utah Raster toolkit and
X# has provisions for displaying a quantized image on an Iris workstation.
X#
X# Location of Utah Raster include directory and library:
XRLEINC = /usr/u/utah/include
XRLELIB = /usr/u/utah/lib/librle.a
X#
X# If you've got an IRIS 4D, use the following two lines.
XCFLAGS = -O -I$(RLEINC) -DIRIS -float
XLIBS = $(RLELIB) -lgl_s -lm
X#
X# Else, use some variation on the following two lines.
X#CFLAGS = -O -I$(RLEINC)
X#LIBS = $(RLELIB) -lm
X
Xall: quantdemo
X
Xquantdemo:	quantdemo.o colorquant.o
X		cc $(CFLAGS) -o quantdemo quantdemo.o colorquant.o $(LIBS)
X
Xclean:
X	/bin/rm -f *.o
END_OF_FILE
if test 711 -ne `wc -c <'Makefile'`; then
    echo shar: \"'Makefile'\" unpacked with wrong size!
fi
# end of 'Makefile'
fi
if test -f 'PORTING' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'PORTING'\"
else
echo shar: Extracting \"'PORTING'\" \(1170 characters\)
sed "s/^X//" >'PORTING' <<'END_OF_FILE'
XThe colorquant() routine has been tested on an Iris 4D70 workstation,
Xa Sun 3/60 workstation, and (to some extent), a Macintosh.
X
XCalls to bzero() may have to be replaced with the appropriate thing on
Xyour system.  (bzero(ptr, len) writes 'len' 0-bytes starting at the location
Xpointed to by ptr.)
X
XAlthough I've tried to avoid integer overflow problems where ever possible,
Xit's likely I've missed a spot where an 'int' should really be a 'long'.
X(On the machine this was developed on, an int == long == 32 bits.)
X
XNote that it's quite easy to optimize this code for a given value for
X'bits'.  In addition, if you assume bits is small and
Xthat the total number of pixels is relatively small, there are several
Xplaces that integer arithmetic may be substituted for floating-point.
XOne such place is the loop in BoxStats -- mean and var need not necessary
Xbe floats.
X
XAs things stand, the maximum number of colors to which an image may
Xbe quantized is 256.  This limit may be overcome by changing rgbmap and
Xcolormap from arrays of characters to arrays of something larger.
X
XSee 'quantdemo.c' for an example of how colorquant() is used.
X
XCraig Kolb 8/25/89
Xkolb@yale.edu
END_OF_FILE
if test 1170 -ne `wc -c <'PORTING'`; then
    echo shar: \"'PORTING'\" unpacked with wrong size!
fi
# end of 'PORTING'
fi
if test -f 'colorquant.3' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'colorquant.3'\"
else
echo shar: Extracting \"'colorquant.3'\" \(2983 characters\)
sed "s/^X//" >'colorquant.3' <<'END_OF_FILE'
X.TH COLORQUANT 3 "August 14, 1989"
X.UC 4
X.SH NAME
Xcolorquant \- variance-based color quantization
X.SH SYNOPSIS
Xint colorquant(red, green, blue, npix, colormap, colors, bits, rgbmap, fast)
X.br
Xunsigned char *red, *green, *blue;
X.br
Xunsigned long npix;
X.br
Xunsigned char *colormap[3];
X.br
Xint colors, bits;
X.br
Xunsigned char *rgbmap;
X.br
Xint fast;
X.SH DESCRIPTION
X.I Colorquant 
Xperforms variance-based color quantization on a given image.
XA representative colormap
Xand a table for performing RGB to colormap index mapping are computed.  The
Xnumber of colors to which the image was quantized (the total number
Xof colormap entries computed) is returned.
XThe arguments to
X.I colorquant 
Xare:
X.TP
X.B red, green, blue
XThe red, green and blue channels of the image.  The ith pixel is represented
Xas the RGB triple (\fBred\fR[i], \fBgreen\fR[i], \fBblue\fR[i]).  These
Xarrays usually contain values which have been 'prequantized' (see below).
X.TP
X.B npix
XThe length, in bytes, of the \fBred\fR, \fBgreen\fR and \fBblue\fR arrays.
XIn other words, the total number of pixels in the image.
X.TP
X.B colormap
XPoints to a pre-allocated, three-channel colormap.  These arrays will be
Xfilled with the colormap values computed by the variance-based color
Xquantization algorithm.  \fBcolormap\fR[0][i], \fBcolormap\fR[1][i], and
X\fBcolormap\fR[2][i] are, respectively, the red, green and blue components
Xof the ith colormap entry.
X.TP
X.B colors
XThe number of pre-allocated colormap entries.  The image will be quantized to
Xat most this many colors.
X.TP
X.B bits
XThe number of significant bits in each entry of the \fBred\fR, \fBgreen\fR and
X\fBblue\fR arrays.  Normally, the red, green and blue arrays contain
Xvalues which have been prequantized from eight to a lower number of
Xsignificant bits.
XFive significant bits usually represents a good tradeoff between image quality
Xand running time.  Anything above six significant bits will likely lead to
Xexcessive paging, as the size of \fBrgbmap\fR and the internal histogram are
Xproportional to (2^\fBbits\fR)^3. 
X.TP
X.B rgbmap
XA pointer to an array of unsigned chars of size (2^\fBbits\fR)^3.
XThis array is used
Xto map from pixels to colormap entries.  The prequantized red, green
Xand blue components of a pixel are used as an index into this array
Xto retrieve the colormap index which should be used to represent the
Xpixel.  The array is indexed as:
X.ce 1
Xcolorindex = \fBrgbmap\fR[(((r << \fBbits\fR) | g) << \fBbits\fR) | b];
Xwhere r, g, and b are the prequantized red, green and blue components of
Xthe pixel in question.
X.TP
X.B fast
XIf non-zero, the construction of rgbmap will be relatively fast.  If
Xzero, \fBrgbmap\fR will be built slowly but more accurately.  In most cases,
Xthe error introduced by the 'fast' approximation is barely noticeable.
X.SH AUTHOR
XCraig Kolb
X.SH REFERENCE
XWan, Wong, and Prusinkiewicz,
X\fIAn Algorithm for Multidimensional Data Clustering,\fR
XTransactions on Mathematical Software, Vol. 14 #2 (June, 1988), pp. 153-162.
END_OF_FILE
if test 2983 -ne `wc -c <'colorquant.3'`; then
    echo shar: \"'colorquant.3'\" unpacked with wrong size!
fi
# end of 'colorquant.3'
fi
if test -f 'colorquant.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'colorquant.c'\"
else
echo shar: Extracting \"'colorquant.c'\" \(16211 characters\)
sed "s/^X//" >'colorquant.c' <<'END_OF_FILE'
X/*
X * This software is copyrighted as noted below.  It may be freely copied,
X * modified, and redistributed, provided that the copyright notice is 
X * preserved on all copies.
X * 
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely "as is".  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the 
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X */
X/*
X * colorquant.c
X *
X * Perform variance-based color quantization on a "full color" image.
X * Author:	Craig Kolb
X *		Department of Mathematics
X *		Yale University
X *		kolb@yale.edu
X * Date:	Tue Aug 22 1989
X * Copyright (C) 1989 Craig E. Kolb
X */
X#include <math.h>
X#include <stdio.h>
X
X/*
X * Maximum number of colormap entries.  To make larger than 2^8, the rgbmap
X * type will have to be changed from unsigned chars to something larger.
X */
X#define MAXCOLORS		256
X/*
X * Value corresponding to full intensity in colormap.  The values placed
X * in the colormap are scaled to be between zero and this number.  Note
X * that anything larger than 255 is going to lead to problems, as the
X * colormap is declared as an unsigned char.
X */
X#define FULLINTENSITY		255
X#define MAX(x,y)	((x) > (y) ? (x) : (y))
X
X/*
X * Readability constants.
X */
X#define REDI		0	
X#define GREENI		1
X#define BLUEI		2	
X#define TRUE		1
X#define FALSE		0
X
Xtypedef struct {
X	float		weightedvar,		/* weighted variance */
X			mean[3];		/* centroid */
X	unsigned long 	weight,			/* # of pixels in box */
X			freq[3][MAXCOLORS];	/* Projected frequencies */
X	int 		low[3], high[3];	/* Box extent */
X} Box;
X
Xunsigned long		*Histogram,		/* image histogram */
X			NPixels;		/* total # of pixels */
Xunsigned int		Colors,			/* Maximum colormap entries */
X			Bits,			/* # significant input bits */
X			ColormaxI;		/* # of colors, 2^Bits */
X
X/*
X * Perform variance-based color quantization on a 24-bit image.
X *
X * Input consists of:
X *	red, green, blue	Arrays of red, green and blue pixel
X *				intensities stored as unsigned characters.
X *				The color of the ith pixel is given
X *				by red[i] green[i] and blue[i].  0 indicates
X *				zero intensity, 255 full intensity.
X *	pixels			The length of the red, green and blue arrays
X *				in bytes, stored as an unsigned long int.
X *	colormap		Points to the colormap.  The colormap
X *				consists of red, green and blue arrays.
X *				The red/green/blue values of the ith
X *				colormap entry are given respectively by
X *				colormap[0][i], colormap[1][i] and
X *				colormap[2][i].  Each entry is an unsigned char.
X *	colors			The number of colormap entries, stored
X *				as an integer.	
X *	bits			The number of significant bits in each entry
X *				of the red, green and blue arrays. An integer.
X *	rgbmap			An array of unsigned chars of size (2^bits)^3.
X *				This array is used to map from pixels to
X *				colormap entries.  The 'prequantized' red,
X *				green and blue components of a pixel
X *				are used as an index into rgbmap to retrieve
X *				the index which should be used into the colormap
X *				to represent the pixel.  In short:
X *				index = rgbmap[(((r << bits) | g) << bits) | b];
X * 	fast			If non-zero, the rgbmap will be constructed
X *				quickly.  If zero, the rgbmap will be built
X *				much slower, but more accurately.  In most
X *				cases, fast should be non-zero, as the error
X *				introduced by the approximation is usually
X *				small.  'Fast' is stored as an integer.
X *
X * colorquant returns the number of colors to which the image was
X * quantized.
X */
Xint
Xcolorquant(red,green,blue,pixels,colormap,colors,bits,rgbmap,fast)
Xunsigned char *red, *green, *blue;
Xunsigned long pixels;
Xunsigned char *colormap[3];
Xint colors, bits;
Xunsigned char *rgbmap;
Xint fast;
X{
X	Box	*Boxes;				/* Array of color boxes. */
X	int	i,				/* Counter */
X		OutColors,			/* # of entries computed */
X		Colormax;			/* quantized full-intensity */ 
X	float	Cfactor;			/* Conversion factor */
X
X	ColormaxI = 1 << bits;			/* 2 ^ Bits */
X	Colormax = ColormaxI - 1;
X	Bits = bits;
X	NPixels = pixels;
X	Cfactor = (float)FULLINTENSITY / Colormax;
X
X	Histogram = (unsigned long *)calloc(ColormaxI*ColormaxI*ColormaxI, 
X						sizeof(long));
X	Boxes = (Box *)malloc(colors * sizeof(Box));
X	QuantHistogram(red, green, blue, &Boxes[0]);
X	OutColors = CutBoxes(Boxes, colors);
X	/*
X	 * We now know the set of representative colors.  We now
X	 * must fill in the colormap and convert the representatives
X	 * from their 'prequantized' range to 0-FULLINTENSITY.
X	 */
X	for (i = 0; i < OutColors; i++) {
X		colormap[0][i] =
X			(unsigned char)(Boxes[i].mean[REDI] * Cfactor + 0.5);
X		colormap[1][i] =
X			(unsigned char)(Boxes[i].mean[GREENI] * Cfactor + 0.5);
X		colormap[2][i] =
X			(unsigned char)(Boxes[i].mean[BLUEI] * Cfactor + 0.5);
X	}
X
X	ComputeRGBMap(Boxes, OutColors, rgbmap, bits, colormap, fast);
X	free((char *)Histogram);
X	free((char *)Boxes);
X
X	return OutColors;		/* Return # of colormap entries */
X}
X
X/*
X * Compute the histogram of the image as well as the projected frequency
X * arrays for the first world-encompassing box.
X */
XQuantHistogram(r, g, b, box)
Xregister unsigned char *r, *g, *b;
XBox *box;
X{
X	unsigned long *rf, *gf, *bf, i;
X
X	rf = box->freq[0];
X	gf = box->freq[1];
X	bf = box->freq[2];
X
X	/*
X	 * Zero-out the projected frequency arrays of the largest box.
X	 * We compute both the histogram and the proj. frequencies of
X	 * the first box at the same time to save a pass through the
X	 * entire image. 
X	 */
X	bzero(rf, ColormaxI * sizeof(unsigned long));
X	bzero(gf, ColormaxI * sizeof(unsigned long));
X	bzero(bf, ColormaxI * sizeof(unsigned long));
X
X	for (i = 0; i < NPixels; i++) {
X		rf[*r]++;
X		gf[*g]++;
X		bf[*b]++;
X		Histogram[((((*r++)<<Bits)|(*g++))<<Bits)|(*b++)]++;
X	}
X}
X
X/*
X * Iteratively cut the boxes.
X */
XCutBoxes(boxes, colors) 
XBox	*boxes;
Xint	colors;
X{
X	int curbox;
X
X	boxes[0].low[REDI] = boxes[0].low[GREENI] = boxes[0].low[BLUEI] = 0;
X	boxes[0].high[REDI] = boxes[0].high[GREENI] =
X			      boxes[0].high[BLUEI] = ColormaxI;
X	boxes[0].weight = NPixels;
X
X	BoxStats(&boxes[0]);
X
X	for (curbox = 1; curbox < colors; curbox++) {
X		if (CutBox(&boxes[GreatestVariance(boxes, curbox)],
X			   &boxes[curbox]) == FALSE)
X				break;
X	}
X
X	return curbox;
X}
X
X/*
X * Return the number of the box in 'boxes' with the greatest variance.
X * Restrict the search to those boxes with indices between 0 and n-1.
X */
XGreatestVariance(boxes, n)
XBox *boxes;
Xint n;
X{
X	register int i, whichbox;
X	float max;
X
X	max = -1;
X	for (i = 0; i < n; i++) {
X		if (boxes[i].weightedvar > max) {
X			max = boxes[i].weightedvar;
X			whichbox = i;
X		}
X	}
X	return whichbox;
X}
X
X/*
X * Compute mean and weighted variance of the given box.
X */
XBoxStats(box)
Xregister Box *box;
X{
X	register int i, color;
X	unsigned long *freq;
X	float mean, var;
X
X	if(box->weight == 0) {
X		box->weightedvar = 0;
X		return;
X	}
X
X	box->weightedvar = 0.;
X	for (color = 0; color < 3; color++) {
X		var = mean = 0;
X		i = box->low[color];
X		freq = &box->freq[color][i];
X		for (; i < box->high[color]; i++, freq++) {
X			mean += i * *freq;
X			var += i*i* *freq;
X		}
X		box->mean[color] = mean / (float)box->weight;
X		box->weightedvar += var - box->mean[color]*box->mean[color]*
X					(float)box->weight;
X	}
X	box->weightedvar /= NPixels;
X}
X
X/*
X * Cut the given box.  Returns TRUE if the box could be cut, FALSE otherwise.
X */
XCutBox(box, newbox)
XBox *box, *newbox;
X{
X	int i;
X	float totalvar[3];
X	Box newboxes[3][2];
X
X	if (box->weightedvar == 0. || box->weight == 0)
X		/*
X		 * Can't cut this box.
X		 */
X		return FALSE;
X
X	/*
X	 * Find 'optimal' cutpoint along each of the red, green and blue
X	 * axes.  Sum the variances of the two boxes which would result
X	 * by making each cut and store the resultant boxes for 
X	 * (possible) later use.
X	 */
X	for (i = 0; i < 3; i++) {
X		FindCutpoint(box, i, &newboxes[i][0], &newboxes[i][1]);
X		totalvar[i] = newboxes[i][0].weightedvar +
X				newboxes[i][1].weightedvar;
X	}
X
X	/*
X	 * Find which of the three cuts minimized the total variance
X	 * and make that the 'real' cut.
X	 */
X	if (totalvar[REDI] <= totalvar[GREENI] &&
X	    totalvar[REDI] <= totalvar[BLUEI]) {
X		*box = newboxes[REDI][0];
X		*newbox = newboxes[REDI][1];
X	} else if (totalvar[GREENI] <= totalvar[REDI] &&
X		 totalvar[GREENI] <= totalvar[BLUEI]) {
X		*box = newboxes[GREENI][0];
X		*newbox = newboxes[GREENI][1];
X	} else {
X		*box = newboxes[BLUEI][0];
X		*newbox = newboxes[BLUEI][1];
X	}
X
X	return TRUE;
X}
X
X/*
X * Compute the 'optimal' cutpoint for the given box along the axis
X * indicated by 'color'.  Store the boxes which result from the cut
X * in newbox1 and newbox2.
X */
XFindCutpoint(box, color, newbox1, newbox2)
XBox *box, *newbox1, *newbox2;
Xint color;
X{
X	float u, v, max;
X	int i, maxindex, minindex, cutpoint;
X	unsigned long optweight, curweight;
X
X	if (box->low[color] + 1 == box->high[color]) {
X		newbox1->weightedvar = newbox2->weightedvar = HUGE;
X		return;
X	}
X	minindex = (int)((box->low[color] + box->mean[color]) * 0.5);
X	maxindex = (int)((box->mean[color] + box->high[color]) * 0.5);
X
X	cutpoint = minindex;
X	optweight = box->weight;
X
X	curweight = 0.;
X	for (i = box->low[color] ; i < minindex ; i++)
X		curweight += box->freq[color][i];
X	u = 0.;
X	max = -1;
X	for (i = minindex; i <= maxindex ; i++) {
X		curweight += box->freq[color][i];
X		if (curweight == box->weight)
X			break;
X		u += (float)(i * box->freq[color][i]) /
X					(float)box->weight;
X		v = ((float)curweight / (float)(box->weight-curweight)) *
X				(box->mean[color]-u)*(box->mean[color]-u);
X		if (v > max) {
X			max = v;
X			cutpoint = i;
X			optweight = curweight;
X		}
X	}
X	cutpoint++;
X	*newbox1 = *newbox2 = *box;
X	newbox1->weight = optweight;
X	newbox2->weight -= optweight;
X	newbox1->high[color] = cutpoint;
X	newbox2->low[color] = cutpoint;
X	UpdateFrequencies(newbox1, newbox2);
X	BoxStats(newbox1);
X	BoxStats(newbox2);
X}
X
X/*
X * Update projected frequency arrays for two boxes which used to be
X * a single box.
X */
XUpdateFrequencies(box1, box2)
Xregister Box *box1, *box2;
X{
X	register unsigned long myfreq, *h;
X	register int b, g, r;
X	int roff;
X
X	bzero(box1->freq[0], ColormaxI * sizeof(unsigned long));
X	bzero(box1->freq[1], ColormaxI * sizeof(unsigned long));
X	bzero(box1->freq[2], ColormaxI * sizeof(unsigned long)); 
X
X	for (r = box1->low[0]; r < box1->high[0]; r++) {
X		roff = r << Bits;
X		for (g = box1->low[1];g < box1->high[1]; g++) {
X			b = box1->low[2];
X			h = Histogram + (((roff | g) << Bits) | b);
X			for (; b < box1->high[2]; b++) {
X				if ((myfreq = *h++) == 0)
X					continue;
X				box1->freq[0][r] += myfreq;
X				box1->freq[1][g] += myfreq;
X				box1->freq[2][b] += myfreq;
X				box2->freq[0][r] -= myfreq;
X				box2->freq[1][g] -= myfreq;
X				box2->freq[2][b] -= myfreq;
X			}
X		}
X	}
X}
X
X/*
X * Compute RGB to colormap index map.
X */
XComputeRGBMap(boxes, colors, rgbmap, bits, colormap, fast)
XBox *boxes;
Xint colors;
Xunsigned char *rgbmap, *colormap[3];
Xint bits, fast;
X{
X	register int i;
X
X	if (fast) {
X		/*
X		 * The centroid of each box serves as the representative
X		 * for each color in the box.
X		 */
X		for (i = 0; i < colors; i++)
X			SetRGBmap(i, &boxes[i], rgbmap, bits);
X	} else
X		/*
X		 * Find the 'nearest' representative for each
X		 * pixel.
X		 */
X		find_colors(boxes, colors, rgbmap, bits, colormap);
X}
X
X/*
X * Make the centroid of "boxnum" serve as the representative for
X * each color in the box.
X */
XSetRGBmap(boxnum, box, rgbmap, bits)
Xint boxnum;
XBox *box;
Xunsigned char *rgbmap;
Xint bits;
X{
X	register int r, g, b;
X	
X	for (r = box->low[REDI]; r < box->high[REDI]; r++) {
X		for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
X			for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
X				rgbmap[(((r<<bits)|g)<<bits)|b]=(char)boxnum;
X			}
X		}
X	}
X}
X/*
X * Form colormap and NearestColor arrays.
X */
Xfind_colors(boxes, colors, rgbmap, bits, colormap)
Xint colors;
XBox *boxes;
Xunsigned char *rgbmap;
Xint bits;
Xunsigned char *colormap[3];
X{
X	register int i;
X	int num, *neighbors;
X
X	neighbors = (int *)malloc(colors * sizeof(int));
X
X	/*
X	 * Form map of representative (nearest) colors.
X	 */
X	for (i = 0; i < colors; i++) {
X		/*
X		 * Create list of candidate neighbors and
X		 * find closest representative for each
X		 * color in the box.
X		 */
X		num = getneighbors(boxes, i, neighbors, colors, colormap);
X		makenearest(boxes, i, num, neighbors, rgbmap, bits, colormap);
X	}
X	free((char *)neighbors);
X}
X
X/*
X * In order to minimize our search for 'best representative', we form the
X * 'neighbors' array.  This array contains the number of the boxes whose
X * centroids *might* be used as a representative for some color in the
X * current box.  We need only consider those boxes whose centroids are closer
X * to one or more of the current box's corners than is the centroid of the
X * current box. 'Closeness' is measured by Euclidean distance.
X */
Xgetneighbors(boxes, num, neighbors, colors, colormap)
XBox *boxes;
Xint num, colors, *neighbors;
Xunsigned char *colormap[3];
X{
X	register int i, j;
X	Box *bp;
X	float dist, LowR, LowG, LowB, HighR, HighG, HighB, ldiff, hdiff;
X
X	bp = &boxes[num];
X
X	ldiff = bp->low[REDI] - bp->mean[REDI];
X	ldiff *= ldiff;
X	hdiff = bp->high[REDI] - bp->mean[REDI];
X	hdiff *= hdiff;
X	dist = MAX(ldiff, hdiff);
X
X	ldiff = bp->low[GREENI] - bp->mean[GREENI];
X	ldiff *= ldiff;
X	hdiff = bp->high[GREENI] - bp->mean[GREENI];
X	hdiff *= hdiff;
X	dist += MAX(ldiff, hdiff);
X
X	ldiff = bp->low[BLUEI] - bp->mean[BLUEI];
X	ldiff *= ldiff;
X	hdiff = bp->high[BLUEI] - bp->mean[BLUEI];
X	hdiff *= hdiff;
X	dist += MAX(ldiff, hdiff);
X
X#ifdef IRIS
X	dist = fsqrt(dist);
X#else
X	dist = (float)sqrt((double)dist);
X#endif
X
X	/*
X	 * Loop over all colors in the colormap, the ith entry of which
X	 * corresponds to the ith box.
X	 *
X	 * If the centroid of a box is as close to any corner of the
X	 * current box as is the centroid of the current box, add that
X	 * box to the list of "neighbors" of the current box.
X	 */
X	HighR = (float)bp->high[REDI] + dist;
X	HighG = (float)bp->high[GREENI] + dist;
X	HighB = (float)bp->high[BLUEI] + dist;
X	LowR = (float)bp->low[REDI] - dist;
X	LowG = (float)bp->low[GREENI] - dist;
X	LowB = (float)bp->low[BLUEI] - dist;
X	for (i = j = 0, bp = boxes; i < colors; i++, bp++) {
X		if (LowR <= bp->mean[REDI] && HighR >= bp->mean[REDI] &&
X		    LowG <= bp->mean[GREENI] && HighG >= bp->mean[GREENI] &&
X		    LowB <= bp->mean[BLUEI] && HighB >= bp->mean[BLUEI])
X			neighbors[j++] = i;
X	}
X
X	return j;	/* Return the number of neighbors found. */
X}
X
X/*
X * Assign representative colors to every pixel in a given box through
X * the construction of the NearestColor array.  For each color in the
X * given box, we look at the list of neighbors passed to find the
X * one whose centroid is closest to the current color.
X */
Xmakenearest(boxes, boxnum, nneighbors, neighbors, rgbmap, bits, colormap)
XBox *boxes;
Xint boxnum;
Xint nneighbors, *neighbors, bits;
Xunsigned char *rgbmap, *colormap[3];
X{
X	register int n, b, g, r;
X	float rdist, gdist, bdist, dist, mindist;
X	int which, *np;
X	Box *box;
X	extern unsigned long *Histogram;
X
X	box = &boxes[boxnum];
X
X	for (r = box->low[REDI]; r < box->high[REDI]; r++) {
X		for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
X			for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
X/*
X * The following two lines should be commented out if the RGBmap is going
X * to be used for images other than the one given.
X */
X				if (Histogram[(((r<<bits)|g)<<bits)|b] == 0)
X					continue;
X				mindist = HUGE;
X				/*
X				 * Find the colormap entry which is
X				 * closest to the current color.
X				 */
X				np = neighbors;
X				for (n = 0; n < nneighbors; n++, np++) {
X					rdist = r-(int)boxes[*np].mean[REDI];
X					gdist = g-(int)boxes[*np].mean[GREENI];
X					bdist = b-(int)boxes[*np].mean[BLUEI];
X					dist = rdist*rdist + gdist*gdist + bdist*bdist;
X					if (dist < mindist) {
X						mindist = dist;
X						which = *np; 
X					}
X				}
X				/*
X				 * The colormap entry closest to this
X				 * color is used as a representative.
X				 */
X				rgbmap[(((r<<bits)|g)<<bits)|b] = which;
X			}
X		}
X	}
X}
END_OF_FILE
if test 16211 -ne `wc -c <'colorquant.c'`; then
    echo shar: \"'colorquant.c'\" unpacked with wrong size!
fi
# end of 'colorquant.c'
fi
if test -f 'quantdemo.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'quantdemo.c'\"
else
echo shar: Extracting \"'quantdemo.c'\" \(6080 characters\)
sed "s/^X//" >'quantdemo.c' <<'END_OF_FILE'
X/*
X * This software is copyrighted as noted below.  It may be freely copied,
X * modified, and redistributed, provided that the copyright notice is 
X * preserved on all copies.
X * 
X * There is no warranty or other guarantee of fitness for this software,
X * it is provided solely "as is".  Bug reports or fixes may be sent
X * to the author, who may or may not act on them as he desires.
X *
X * You may not include this software in a program or other software product
X * without supplying the source, or without informing the end-user that the 
X * source is available for no extra charge.
X *
X * If you modify this software, you should include a notice giving the
X * name of the person performing the modification, the date of modification,
X * and the reason for such modification.
X */
X/*
X * quantdemo.c
X *
X * Demo program for variance-based color quantization routine.
X * Reads a Utah Raster format file and quantizes the image.
X * If compiled with IRIS is defined, code is included to display the
X * quantized image on an Iris 4D workstation.
X *
X * Author:	Craig Kolb
X *		Dept. of Mathematics
X *		Yale University
X *		kolb@yale.edu
X * Date: Wed Feb 8 1989
X * Copyright (C) 1989, Craig E. Kolb
X */
X#include <stdio.h>
X#include <svfb_global.h>
X#include <sys/types.h>
X#include <sys/times.h>
X#include <sys/param.h>
X
X/*
X * DISCLAIMER:
X * This code is sloppy.
X */
X#define DEFAULT_BITS		5
X#define DEFAULT_ENTRIES		256
X#define DEFAULT_FAST		1
X
X#ifdef IRIS
X#include <device.h>
X#define MAP_OFFSET	512	/* Offset into Iris colormap */
Xunsigned short *scanline;
X#endif
X
X#define Quantize(x)	(x >> shift)
X
Xunsigned char *RGBmap, *Pixels[3];	/* RGB-->index map & original image */
Xchar *Progname;				/* argv[0] */
Xint Xres, Yres, Bits;			/* Image resolution, # of sig. bits */
Xint Verbose;				/* Verbosity flag. */
X
Xmain(argc, argv)
Xint argc;
Xchar **argv;
X{
X	unsigned char *rp, *gp, *bp, *inline[3], *colormap[3];
X	unsigned ByteCount;
X	int x, y, colors, Entries, shift, Fast;
X	struct tms time;
X
X
X#ifdef IRIS
X	short val;
X#endif
X	/*
X	 * Set default values.
X	 */
X	Bits = DEFAULT_BITS;
X	Entries = DEFAULT_ENTRIES;
X	Fast = DEFAULT_FAST;
X
X	Progname = argv[0];
X
X	argv++;
X	argc--;
X
X	while (argc) {
X		if (argv[0][0] != '-')
X			break;
X		switch(argv[0][1]) {
X			case 'b':
X				Bits = atoi(argv[1]);
X				argc--;
X				argv++;
X				break;
X			case 'c':
X				Entries = atoi(argv[1]);
X				argc--;
X				argv++;
X				break;
X			case 'f':
X				Fast = !Fast;
X				break;
X			case 'h':
X				usage();
X				exit(3);
X			case 'v':
X				Verbose = !Verbose;
X				break;
X			default:
X				fprintf(stderr,"Unknown option: %s\n",argv[0]);
X				usage();
X				exit(1);
X		}
X		argv++;
X		argc--;
X	}
X
X	if (argc) {
X		sv_globals.svfb_fd = fopen(argv[0], "r");
X		if (sv_globals.svfb_fd == (FILE *)NULL) {
X			fprintf(stderr,"Cannot open %s.\n",argv[0]);
X			exit(0);
X		}
X	} else
X		sv_globals.svfb_fd = stdin;
X
X	/*
X	 * Read Utah Raster header.
X	 */
X	if (rle_get_setup(&sv_globals) < 0) {
X		fprintf(stderr,"Error reading RLE header.\n");
X		exit(1);
X	}
X
X	/*
X	 * Ignore ALPHA channel, if any.
X	 */
X	SV_CLR_BIT(sv_globals, SV_ALPHA);
X	Yres = sv_globals.sv_ymax - sv_globals.sv_ymin + 1;
X	Xres = sv_globals.sv_xmax - sv_globals.sv_xmin + 1;
X	sv_globals.sv_xmax -= sv_globals.sv_xmin;
X	sv_globals.sv_xmin = 0;
X
X	/*
X	 * ByteCount gives the size off the RGBmap, and is equal to
X	 * (2^Bits)^3
X	 */
X	ByteCount = 1 << Bits;
X	ByteCount *= ByteCount * ByteCount;
X	RGBmap = (unsigned char *)malloc(ByteCount*sizeof(unsigned char));
X
X	Pixels[0]= (unsigned char *)malloc(Xres * Yres * sizeof(unsigned char));
X	Pixels[1]= (unsigned char *)malloc(Xres * Yres * sizeof(unsigned char));
X	Pixels[2]= (unsigned char *)malloc(Xres * Yres * sizeof(unsigned char));
X
X	rp = Pixels[0];
X	gp = Pixels[1];
X	bp = Pixels[2];
X
X	inline[0] = (unsigned char *)malloc(Xres * sizeof(unsigned char));
X	inline[1] = (unsigned char *)malloc(Xres * sizeof(unsigned char));
X	inline[2] = (unsigned char *)malloc(Xres * sizeof(unsigned char));
X
X	shift = 8 - Bits;
X	for (y = 0; y < Yres; y++) {
X		rle_getrow(&sv_globals, inline);
X		for (x = 0; x < Xres; x++) {
X			/*
X			 * Quantize input values.
X			 */
X			*rp++ = Quantize(inline[0][x]);
X			*gp++ = Quantize(inline[1][x]);
X			*bp++ = Quantize(inline[2][x]);
X		}
X	}
X
X	free((char *)inline[0]);
X	free((char *)inline[1]);
X	free((char *)inline[2]);
X
X	colormap[0] = (unsigned char *)malloc(Entries * sizeof(unsigned char));
X	colormap[1] = (unsigned char *)malloc(Entries * sizeof(unsigned char));
X	colormap[2] = (unsigned char *)malloc(Entries * sizeof(unsigned char));
X
X	if (Verbose) {
X		times(&time);
X		printf("Preprocessing time: %.2fu  %.2fs\n",
X			(float)time.tms_utime / (float)HZ,
X			(float)time.tms_stime / (float)HZ);
X	}
X
X	colors = colorquant(Pixels[0], Pixels[1], Pixels[2],
X				(unsigned long)(Xres*Yres),
X				colormap, Entries, Bits, RGBmap, Fast);
X	if (Verbose) {
X		times(&time);
X		printf("Total time: %.2fu  %.2fs\n",
X			(float)time.tms_utime / (float)HZ,
X			(float)time.tms_stime / (float)HZ);
X		fprintf(stderr,"Quantized to %d colors.\n",colors);
X	}
X
X#ifdef IRIS
X	prefsize(Xres, Yres);
X	winopen("QuantDemo");
X	unqdevice(INPUTCHANGE);
X	qdevice(REDRAW);
X	qdevice(ESCKEY);
X	qdevice(WINQUIT);
X	/*
X	 * Load hardware colormap with computed values.
X	 */
X	for (x = 0; x < colors; x++)
X		mapcolor(x+MAP_OFFSET, (short)(colormap[0][x]),
X				(short)(colormap[1][x]),
X				(short)(colormap[2][x]));
X	scanline = (unsigned short *)malloc(Xres * sizeof(unsigned short));
X
X	drawimage();
X
X	for (;;) {
X		switch(qread(&val)) {
X			case WINQUIT:
X			case ESCKEY:
X				gexit(0);
X				exit(0);
X				break;
X			case REDRAW:
X				reshapeviewport();
X				drawimage();
X				break;
X		}
X	}
X#else
X	exit(0);
X#endif
X}
X
X#ifdef IRIS
X/*
X * Display quantized image.
X */
Xdrawimage()
X{
X	int x, y;
X	unsigned char *rp, *gp, *bp;
X
X	rp = Pixels[0];
X	gp = Pixels[1];
X	bp = Pixels[2];
X
X	for (y = 0; y < Yres; y++) {
X		for (x = 0; x < Xres; x++, rp++, gp++, bp++)
X			scanline[x] = MAP_OFFSET + (unsigned short)
X				RGBmap[(((*rp<<Bits)|*gp)<<Bits)|*bp];
X		cmov2i(0, y);
X		writepixels(Xres, scanline);
X	}
X}
X#endif
Xusage()
X{
X	fprintf(stderr,"usage: quantdemo [-c colors] [-b bits] [-f] [filename]\n");
X}
END_OF_FILE
if test 6080 -ne `wc -c <'quantdemo.c'`; then
    echo shar: \"'quantdemo.c'\" unpacked with wrong size!
fi
# end of 'quantdemo.c'
fi
echo shar: End of archive 1 \(of 1\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have the archive.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0