/*
 * This software is copyrighted as noted below.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is 
 * preserved on all copies.
 * 
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the author, who may or may not act on them as he desires.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the 
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 */
/*
 * colorquant.c
 *
 * Perform variance-based color quantization on a "full color" image.
 * Author:	Craig Kolb
 *		Department of Mathematics
 *		Yale University
 *		kolb@yale.edu
 * Date:	Tue Aug 22 1989
 * Copyright (C) 1989 Craig E. Kolb
 * $Id: colorquant.c,v 3.0.1.4 1992/04/30 14:06:48 spencer Exp $
 *
 * $Log: /fw4/app/common/corec/colorquant.cpp $
 * 
 * 7     5/8/00 11:11 AM Mbendiksen
 * 
 * 6     5/4/00 11:18 AM Mbendiksen
 * 
 * 5     5/3/00 12:20 PM Mbendiksen
 * 
 * 4     3/20/00 9:50a Sjohnson
 * 
 * 3     1/31/00 4:41p Sjohnson
 * 
 * 2     12/14/99 3:50p Sjohnson
 * 
 * 1     11/30/99 1:31 PM Mbendiksen
 * Adding subproject 'app' to '$/FW4'
 * 
 * 14    2/27/98 10:05 AM mbendiksen
 * 
 * 13    2/22/98 7:22 PM sjohnson
 * 
 * 12    2/09/98 6:19p jeffd
 * 
 * 11    2/09/98 6:07p jeffa
 * 
 * 10    2/9/98 5:23 PM mbendiksen
 * 
 * 9     2/5/98 9:14 AM sjohnson
 * 
 * 8     11/20/97 9:24a sjohnson
 * 
 * 7     11/19/97 5:37 PM sjohnson
 * 
 * 6     10/07/97 1:35p jeffa
 * 
 * 5     10/7/97 12:33 PM mbendiksen
 * 
 * 2     10/01/97 9:01a sjohnson
 * Revision 3.0.1.4  1992/04/30  14:06:48  spencer
 * patch3: lint fixes.
 *
 * Revision 3.0.1.3  1991/02/06  16:46:47  spencer
 * patch3: Fix round-off error in variance calculation that sometimes
 * patch3: caused multiple boxes with the same center to be cut.
 * patch3: Change 'fast' argument to 'flags' argument.  Currently:
 * patch3: 	CQ_FAST: same as old fast argument
 * patch3: 	CQ_QUANTIZE: data is not prequantized to bits
 * patch3: 	CQ_NO_RGBMAP: don't build rgbmap.
 *
 * Revision 3.0.1.2  90/11/29  15:18:04  spencer
 * Remove a typo.
 * 
 * 
 * Revision 3.0.1.1  90/11/19  16:59:48  spencer
 * Use inv_cmap instead of find_colors.
 * Changes to process multiple files into one colormap (accum_hist arg).
 * Delete 'otherimages' argument -- unnecessary with faster inv_cmap code.
 * 
 * 
 * Revision 3.0  90/08/03  15:20:11  spencer
 * Establish version 3.0 base.
 * 
 * Revision 1.6  90/07/29  08:06:06  spencer
 * If HUGE isn't defined, make it HUGE_VAL (for ansi).
 * 
 * Revision 1.5  90/07/26  17:25:48  rgb
 * Added a parameter to colorquant for rgbmap construction.
 * 
 * Revision 1.4  90/07/13  14:53:31  spencer
 * Get rid of ARB_ARG sh*t.
 * Change a couple of vars to double so that HUGE won't cause problems.
 * 
 * Revision 1.3  90/06/28  21:42:56  spencer
 * Declare internal functions properly.
 * Delete unused global variable.
 * 
 * Revision 1.2  90/06/28  13:18:56  spencer
 * Make internal functions static.
 * Build entire RGB cube, not just those colors used by this image.  This
 * is so dithering will work.
 * 
 * Revision 1.1  90/06/18  20:45:17  spencer
 * Initial revision
 * 
 * Revision 1.3  89/12/03  18:27:16  craig
 * Removed bogus integer casts in distance calculation in makenearest().
 * 
 * Revision 1.2  89/12/03  18:13:12  craig
 * FindCutpoint now returns FALSE if the given box cannot be cut.  This
 * to avoid overflow problems in CutBox.
 * "whichbox" in GreatestVariance() is now initialized to 0.
 * 
 */
static char rcsid[] = "$Header: /fw4/app/common/corec/colorquant.cpp 7     5/8/00 11:11 AM Mbendiksen $";

#include "MbStdHeaders.h"	// this should come FIRST in every file.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
//#include "rle_config.h"
	#include "colorquant.h"

/* Ansi uses HUGE_VAL. */
#ifndef HUGE
#ifdef HUGE_VAL
#define HUGE HUGE_VAL
#else
#define HUGE MAXFLOAT
#endif
#endif

#define slow_color_map



/* 
 * The colorquant() routine has been tested on an Iris 4D70 workstation,
 * a Sun 3/60 workstation, and (to some extent), a Macintosh.
 * 
 * Calls to bzero() may have to be replaced with the appropriate thing on
 * your system.  (bzero(ptr, len) writes 'len' 0-bytes starting at the location
 * pointed to by ptr.)
 * 
 * Although I've tried to avoid integer overflow problems where ever possible,
 * it's likely I've missed a spot where an 'int' should really be a 'long'.
 * (On the machine this was developed on, an int == long == 32 bits.)
 * 
 * Note that it's quite easy to optimize this code for a given value for
 * 'bits'.  In addition, if you assume bits is small and
 * that the total number of pixels is relatively small, there are several
 * places that integer arithmetic may be substituted for floating-point.
 * One such place is the loop in BoxStats -- mean and var need not necessary
 * be floats.
 * 
 * As things stand, the maximum number of colors to which an image may
 * be quantized is 256.  This limit may be overcome by changing rgbmap and
 * colormap from pointers to characters to pointers to something larger.
 */

/*
 * Maximum number of colormap entries.  To make larger than 2^8, the rgbmap
 * type will have to be changed from unsigned chars to something larger.
 */
#define MAXCOLORS		256
/*
 * Value corrersponding to full intensity in colormap.  The values placed
 * in the colormap are scaled to be between zero and this number.  Note
 * that anything larger than 255 is going to lead to problems, as the
 * colormap is declared as an unsigned char.
 */
#define FULLINTENSITY		255
#ifndef MAX
//This is defined elsewhere on the Mac
#define MAX(x,y)	((x) > (y) ? (x) : (y))
#endif
/*
 * Readability constants.
 */
#define REDI		0	
#define GREENI		1
#define BLUEI		2	

#undef TRUE
#undef FALSE
#define TRUE		1
#define FALSE		0

typedef struct {
	double		weightedvar;		/* weighted variance */
	float		mean[3];			/* centroid */
	unsigned long 	weight,			/* # of pixels in box */
				freq[3][MAXCOLORS];	/* Projected frequencies */
	int 		low[3], high[3];	/* Box extent */
} Box;

unsigned long	*gHistogram = NULL;	/* image histogram */

#ifdef MACINTOSH
Handle			ghHistogramMac = NULL;	// handle to above histogram if alloc'd in temp memory
#endif

unsigned long	gHistogramSize;
long			gNumPixels,			/* # of pixels in an image*/
				gSumPixels;			/* total # of pixels */
static unsigned int	gNumBits,		/* # significant input bits */
				gColormaxI;			/* # of colors, 2^Bits */
static Box		*gBoxes = NULL;		/* Array of color boxes. */


static void
QuantHistogram(unsigned char *r, unsigned char *g, unsigned char *b, 
Box *box, int quantize);
static void BoxStats(Box *box);
static void UpdateFrequencies(Box *box1, Box *box2);
static void SetRGBmap(int boxnum, Box *box, unsigned char *rgbmap, int bits);
#ifdef slow_color_map
static void find_colors(Box *boxes, int colors, unsigned char *rgbmap, int bits,
unsigned char *colormap[3], int otherimages);
static int getneighbors(Box *boxes, int num, int *neighbors,  int colors, 
	unsigned char *colormap[3]);
static void makenearest(Box *boxes, int boxnum, int  nneighbors, int *neighbors, unsigned char *rgbmap, int bits, unsigned char *colormap[3],
	int otherimages);
#else
extern void inv_cmap();
#endif
static int
CutBoxes(Box	*boxes, int	colors);
static int CutBox(Box *box, Box *newbox);
static int GreatestVariance(Box *boxes, int n);
static int FindCutpoint(Box *box, int color, Box *newbox1, Box *newbox2);
/*
 * Perform variance-based color quantization on a 24-bit image.
 *
 * Input consists of:
 *	red, green, blue	Arrays of red, green and blue pixel
 *				intensities stored as unsigned characters.
 *				The color of the ith pixel is given
 *				by red[i] green[i] and blue[i].  0 indicates
 *				zero intensity, 255 full intensity.
 *	pixels			The length of the red, green and blue arrays
 *				in bytes, stored as an unsigned long int.
 *	colormap		Points to the colormap.  The colormap
 *				consists of red, green and blue arrays.
 *				The red/green/blue values of the ith
 *				colormap entry are given respectively by
 *				colormap[0][i], colormap[1][i] and
 *				colormap[2][i].  Each entry is an unsigned char.
 *	colors			The number of colormap entries, stored
 *				as an integer.	
 *	bits			The number of significant bits in each entry
 *				of the red, greena and blue arrays. An integer.
 *	rgbmap			An array of unsigned chars of size (2^bits)^3.
 *				This array is used to map from pixels to
 *				colormap entries.  The 'prequantized' red,
 *				green and blue components of a pixel
 *				are used as an index into rgbmap to retrieve
 *				the index which should be used into the colormap
 *				to represent the pixel.  In short:
 *				index = rgbmap[(((r << bits) | g) << bits) | b];
 * 	flags			A set of bit-flags:
 * 		CQ_FAST		If set, the rgbmap will be constructed
 *				quickly.  If not set, the rgbmap will
 *				be built more slowly, but more
 *				accurately.  In most cases, fast
 *				should be set, as the error introduced
 *				by the approximation is usually small.
 * 		CQ_QUANTIZE:	If set, the data in red, green, and
 * 				blue is not pre-quantized, and will be
 * 				quantized "on the fly".  This slows
 * 				the routine slightly.  If not set, the
 * 				data has been prequantized to 'bits'
 * 				significant bits.
 * 		
 *	accum_hist		If non-zero the histogram will accumulate and 
 *				reflect pixels from multiple images.
 *				when 1, the histogram will be initialized and
 *				summed, but not thrown away OR processed. when 
 *				2 the image RGB will be added to it.  When 3 
 *				Boxes are cut and a colormap and rgbmap
 *				are be returned, Histogram is freed too.
 *				When zero, all code is executed as per normal.
 *
 * colorquant returns the number of colors to which the image was
 * quantized.
 */

/*
 * Compute RGB to colormap index map.
 */
static void
ComputeRGBMap(Box *boxes,
int colors,
unsigned char *rgbmap, int bits, unsigned char *colormap[3],
int fast)
{
	register int i;

	if (fast) {
		/*
		 * The centroid of each box serves as the representative
		 * for each color in the box.
		 */
		for (i = 0; i < colors; i++)
			SetRGBmap(i, &boxes[i], rgbmap, bits);
	} else
		/*
		 * Find the 'nearest' representative for each
		 * pixel.
		 */
#ifdef slow_color_map
		find_colors(boxes, colors, rgbmap, bits, colormap, 1);
#else
		inv_cmap(colors, colormap, bits, gHistogram, rgbmap);
#endif
}


void colorquantCleanup(void) 
{
	if (gHistogram) {
#ifdef MACINTOSH
		if (ghHistogramMac) {
			::HUnlock(ghHistogramMac);
			::DisposeHandle(ghHistogramMac);
			ghHistogramMac = NULL;
			gHistogram = NULL;
		}
		else {
			mb_raw_free((char *)gHistogram);
			gHistogram = NULL;
		}
#else
		mb_raw_free((char *)gHistogram);
		gHistogram = NULL;
#endif
	}
	if (gBoxes) {
		mb_raw_free((char *)gBoxes);
		gBoxes = NULL;
	}
}

int
colorquant(unsigned char *red,unsigned char *green,
		   unsigned char *blue,unsigned long pixels,unsigned char *colormap[3],
		   int colors,int bits,unsigned char *rgbmap,int flags,int accum_hist)
{
    int	i,			/* Counter */
    OutColors,			/* # of entries computed */
    Colormax;			/* quantized full-intensity */ 
    float	Cfactor;	/* Conversion factor */

    if (accum_hist < 0 || accum_hist > PROCESS_HIST)
	fprintf(stderr, "colorquant: bad value for accum_hist\n");

    gColormaxI = 1 << bits;	/* 2 ^ gNumBits */
    Colormax = gColormaxI - 1;
    gNumBits = bits;
    gNumPixels = pixels;
    Cfactor = (float)FULLINTENSITY / Colormax;

    if (!accum_hist || accum_hist == INIT_HIST) {
    	try {	// 980209 added proper failure handling. (mmb)
			gHistogramSize = gColormaxI*gColormaxI*gColormaxI;
			MB_ASSERTCRASH(!gBoxes, ("gBoxes should be NULL"));
			MB_ASSERTCRASH(!gHistogram, ("gHistogram should be NULL"));
			gBoxes = NULL;
			gHistogram = NULL;
#ifdef MACINTOSH	// 980226 (mmb)
			MB_ASSERTCRASH(!ghHistogramMac, ("ghHistogramMac should be NULL"));
			ghHistogramMac = NULL;
			try {
				MB_THROW_NULL(gHistogram = (unsigned long *)mb_raw_malloc(gHistogramSize * sizeof(long)));
			}
			catch (...) {	// failed to alloc in our heap; try temp memory
				OSErr osErr;
				ghHistogramMac = ::TempNewHandle(gHistogramSize * sizeof(long), &osErr);
				if (ghHistogramMac == NULL || osErr != noErr)
					throw EOutOfMem();	// give up!
				::HLock(ghHistogramMac);
				gHistogram = (unsigned long *)(*ghHistogramMac);
				memset((void *)(gHistogram), 0, gHistogramSize * sizeof(long));
			}
#else
			MB_THROW_NULL(gHistogram = (unsigned long *)mb_raw_malloc(gHistogramSize * sizeof(long)));
#endif
			MB_THROW_NULL(gBoxes = (Box *)mb_raw_malloc(colors * sizeof(Box)));
			/*
			 * Zero-out the projected frequency arrays of the largest box.
			 */
			//bzero(gBoxes->freq[0], gColormaxI * sizeof(unsigned long));
			//bzero(gBoxes->freq[1], gColormaxI * sizeof(unsigned long));
			//bzero(gBoxes->freq[2], gColormaxI * sizeof(unsigned long));
			gSumPixels = 0;
    	}
    	catch (...) {
			if (gHistogram) {
#ifdef MACINTOSH
				if (ghHistogramMac) {
					::HUnlock(ghHistogramMac);
					::DisposeHandle(ghHistogramMac);
					ghHistogramMac = NULL;
					gHistogram = NULL;
				}
				else {
					mb_raw_free((char *)gHistogram);
					gHistogram = NULL;
				}
#else
				mb_raw_free((char *)gHistogram);
				gHistogram = NULL;
#endif
			}
			if (gBoxes) {
				mb_raw_free((char *)gBoxes);
				gBoxes = NULL;
			}
			throw;
    	}
    }

    gSumPixels += gNumPixels;

    if ( accum_hist != PROCESS_HIST ) 
	QuantHistogram(red, green, blue, &gBoxes[0], flags&CQ_QUANTIZE);

    if ( !accum_hist || accum_hist == PROCESS_HIST) {
		OutColors = CutBoxes(gBoxes, colors);

		/*
		 * We now know the set of representative colors.  We now
		 * must fill in the colormap and convert the representatives
		 * from their 'prequantized' range to 0-FULLINTENSITY.
		 */
		for (i = 0; i < OutColors; i++) {
			colormap[0][i] =
			(unsigned char)(gBoxes[i].mean[REDI] * Cfactor + 0.5);
			colormap[1][i] =
			(unsigned char)(gBoxes[i].mean[GREENI] * Cfactor + 0.5);
			colormap[2][i] =
			(unsigned char)(gBoxes[i].mean[BLUEI] * Cfactor + 0.5);
		}

		if ( !(flags & CQ_NO_RGBMAP) )
			ComputeRGBMap(gBoxes, OutColors, rgbmap, bits, colormap,
				  flags&CQ_FAST);

		MB_ASSERTCRASH(gHistogram, ("gHistogram should be non-NULL"));
		MB_ASSERTCRASH(gBoxes, ("gBoxes should be non-NULL"));
		if (gHistogram) {
#ifdef MACINTOSH
			if (ghHistogramMac) {
				::HUnlock(ghHistogramMac);
				::DisposeHandle(ghHistogramMac);
				ghHistogramMac = NULL;
				gHistogram = NULL;
			}
			else {
				mb_raw_free((char *)gHistogram);
				gHistogram = NULL;
			}
#else
			mb_raw_free((char *)gHistogram);
			gHistogram = NULL;
#endif
		}
		if (gBoxes) {
			mb_raw_free((char *)gBoxes);
			gBoxes = NULL;
		}
		return OutColors;	/* Return # of colormap entries */
    }

    return 0;
}

/*
 * Compute the histogram of the image as well as the projected frequency
 * arrays for the first world-encompassing box.
 */
static void
QuantHistogram(unsigned char *r, unsigned char *g, unsigned char *b, 
Box *box,
int quantize)
{
	register unsigned long *rf, *gf, *bf;
	long i;

	rf = box->freq[0];
	gf = box->freq[1];
	bf = box->freq[2];

	/*
	 * We compute both the histogram and the proj. frequencies of
	 * the first box at the same time to save a pass through the
	 * entire image. 
	 */
	
	if ( !quantize )
	    for (i = 0; i < gNumPixels; i++) {
		rf[*r]++;
		gf[*g]++;
		bf[*b]++;
		gHistogram[((((*r++)<<gNumBits)|(*g++))<<gNumBits)|(*b++)]++;
	    }
	else
	{
	    unsigned char rr, gg, bb;
	    for (i = 0; i < gNumPixels; i++) {
		rr = (*r++)>>(8-gNumBits);
		gg = (*g++)>>(8-gNumBits);
		bb = (*b++)>>(8-gNumBits);
		rf[rr]++;
		gf[gg]++;
		bf[bb]++;
		gHistogram[(((rr<<gNumBits)|gg)<<gNumBits)|bb]++;
	    }
	}
	    
}

/*
 * Interatively cut the boxes.
 */
static int
CutBoxes(Box	*boxes, int	colors)
{
	int curbox;

	boxes[0].low[REDI] = boxes[0].low[GREENI] = boxes[0].low[BLUEI] = 0;
	boxes[0].high[REDI] = boxes[0].high[GREENI] =
			      boxes[0].high[BLUEI] = gColormaxI;
	boxes[0].weight = gSumPixels;

	BoxStats(&boxes[0]);

	for (curbox = 1; curbox < colors; curbox++) {
		if (CutBox(&boxes[GreatestVariance(boxes, curbox)],
			   &boxes[curbox]) == FALSE)
				break;
	}

	return curbox;
}

/*
 * Return the number of the box in 'boxes' with the greatest variance.
 * Restrict the search to those boxes with indices between 0 and n-1.
 */
static int
GreatestVariance(Box *boxes, int n)
{
	register int i, whichbox = 0;
	double max;

	max = -1;
	for (i = 0; i < n; i++) {
		if (boxes[i].weightedvar > max) {
			max = boxes[i].weightedvar;
			whichbox = i;
		}
	}
	return whichbox;
}

/*
 * Compute mean and weighted variance of the given box.
 */
static void
BoxStats(Box *box)
{
	register int i, color;
	unsigned long *freq;
	double mean, var;

	if(box->weight == 0) {
		box->weightedvar = 0;
		return;
	}

	box->weightedvar = 0.;
	for (color = 0; color < 3; color++) {
		var = mean = 0;
		i = box->low[color];
		freq = &box->freq[color][i];
		for (; i < box->high[color]; i++, freq++) {
			mean += i * *freq;
			var += i*i* *freq;
		}
		box->mean[color] = (float) (mean / (double)box->weight);
		box->weightedvar += var - box->mean[color]*box->mean[color]*
					(double)box->weight;
	}
	box->weightedvar /= gSumPixels;
}

/*
 * Cut the given box.  Returns TRUE if the box could be cut, FALSE otherwise.
 */
static int
CutBox(Box *box, Box *newbox)
{
	int i;
	double totalvar[3];
	Box newboxes[3][2];

	if (box->weightedvar == 0. || box->weight == 0)
		/*
		 * Can't cut this box.
		 */
		return FALSE;

	/*
	 * Find 'optimal' cutpoint along each of the red, green and blue
	 * axes.  Sum the variances of the two boxes which would result
	 * by making each cut and store the resultant boxes for 
	 * (possible) later use.
	 */
	for (i = 0; i < 3; i++) {
		if (FindCutpoint(box, i, &newboxes[i][0], &newboxes[i][1]))
			totalvar[i] = newboxes[i][0].weightedvar +
				newboxes[i][1].weightedvar;
		else
			totalvar[i] = HUGE;
	}

	/*
	 * Find which of the three cuts minimized the total variance
	 * and make that the 'real' cut.
	 */
	if (totalvar[REDI] <= totalvar[GREENI] &&
	    totalvar[REDI] <= totalvar[BLUEI]) {
		*box = newboxes[REDI][0];
		*newbox = newboxes[REDI][1];
	} else if (totalvar[GREENI] <= totalvar[REDI] &&
		 totalvar[GREENI] <= totalvar[BLUEI]) {
		*box = newboxes[GREENI][0];
		*newbox = newboxes[GREENI][1];
	} else {
		*box = newboxes[BLUEI][0];
		*newbox = newboxes[BLUEI][1];
	}

	return TRUE;
}

/*
 * Compute the 'optimal' cutpoint for the given box along the axis
 * indcated by 'color'.  Store the boxes which result from the cut
 * in newbox1 and newbox2.
 */
static int FindCutpoint(Box *box, int color, Box *newbox1, Box *newbox2)
{
	float u, v, max;
	int i, maxindex, minindex, cutpoint;
	unsigned long optweight, curweight;

	if (box->low[color] + 1 == box->high[color])
		return FALSE;	/* Cannot be cut. */
	minindex = (int)((box->low[color] + box->mean[color]) * 0.5);
	maxindex = (int)((box->mean[color] + box->high[color]) * 0.5);

	cutpoint = minindex;
	optweight = box->weight;

	curweight = 0;
	for (i = box->low[color] ; i < minindex ; i++)
		curweight += box->freq[color][i];
	u = 0.;
	max = -1;
	for (i = minindex; i <= maxindex ; i++) {
		curweight += box->freq[color][i];
		if (curweight == box->weight)
			break;
		u += (float)(i * box->freq[color][i]) /
					(float)box->weight;
		v = ((float)curweight / (float)(box->weight-curweight)) *
				(box->mean[color]-u)*(box->mean[color]-u);
		if (v > max) {
			max = v;
			cutpoint = i;
			optweight = curweight;
		}
	}
	cutpoint++;
	*newbox1 = *newbox2 = *box;
	newbox1->weight = optweight;
	newbox2->weight -= optweight;
	newbox1->high[color] = CLIP((long)box->low[color], (long)cutpoint, (long)box->high[color]);
	newbox2->low[color] = CLIP((long)box->low[color], (long)cutpoint, (long)box->high[color]);

	UpdateFrequencies(newbox1, newbox2);
	BoxStats(newbox1);
	BoxStats(newbox2);

	return TRUE;	/* Found cutpoint. */
}

/*
 * Update projected frequency arrays for two boxes which used to be
 * a single box.
 */
static void UpdateFrequencies(Box *box1, Box *box2)
{
	register unsigned long myfreq, *h;
	register int b, g, r;
	int roff;

	{
		unsigned int i;
		for (i=0; i<gColormaxI; i++) {
			box1->freq[0][i] = 0;
			box1->freq[1][i] = 0;
			box1->freq[2][i] = 0;
		}
	}
	
	for (r = box1->low[0]; r < box1->high[0]; r++) {
		roff = r << gNumBits;
		for (g = box1->low[1];g < box1->high[1]; g++) {
			b = box1->low[2];

			long offset = (((roff | g) << gNumBits) | b);
			if (offset >= gHistogramSize) // should never happen but safety first!
				continue;
			h = gHistogram + offset;
			for (; b < box1->high[2]; b++) {
				if ((myfreq = *h++) == 0)
					continue;
				box1->freq[0][r] += myfreq;
				box1->freq[1][g] += myfreq;
				box1->freq[2][b] += myfreq;
				box2->freq[0][r] -= myfreq;
				box2->freq[1][g] -= myfreq;
				box2->freq[2][b] -= myfreq;
			}
		}
	}
}

/*
 * Make the centroid of "boxnum" serve as the representative for
 * each color in the box.
 */
static void SetRGBmap(int boxnum, Box *box, unsigned char *rgbmap, int bits)
{
	register int r, g, b;
	
	for (r = box->low[REDI]; r < box->high[REDI]; r++) {
		for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
			for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
				rgbmap[(((r<<bits)|g)<<bits)|b]=(char)boxnum;
			}
		}
	}
}

#ifdef slow_color_map
/*
 * Form colormap and NearestColor arrays.
 */
static void find_colors(Box *boxes, int colors, unsigned char *rgbmap, int bits,
unsigned char *colormap[3], int otherimages)
{
	register int i;
	int num, *neighbors=NULL;

	try {	// 980209 added proper failure handling. (mmb)
		MB_THROW_NULL(neighbors = (int *)mb_raw_malloc(colors * sizeof(int)));
	
		/*
		 * Form map of representative (nearest) colors.
		 */
		for (i = 0; i < colors; i++) {
			/*
			 * Create list of candidate neighbors and
			 * find closest representative for each
			 * color in the box.
			 */
			num = getneighbors(boxes, i, neighbors, colors, colormap);
			makenearest(boxes, i, num, neighbors, rgbmap, bits, colormap, 
				    otherimages);
		}
		if (neighbors) {
			mb_raw_free((char *)neighbors);
			neighbors = NULL;
		}
	}
	catch (...) {
		if (neighbors) {
			mb_raw_free((char *)neighbors);
			neighbors = NULL;
		}
		throw;
	}
}

/*
 * In order to minimize our search for 'best representative', we form the
 * 'neighbors' array.  This array contains the number of the boxes whose
 * centroids *might* be used as a representative for some color in the
 * current box.  We need only consider those boxes whose centroids are closer
 * to one or more of the current box's corners than is the centroid of the
 * current box. 'Closeness' is measured by Euclidean distance.
 */
static int getneighbors(Box *boxes, int num, int *neighbors,  int colors, 
	unsigned char *colormap[3])
{
	register int i, j;
	Box *bp;
	float dist, LowR, LowG, LowB, HighR, HighG, HighB, ldiff, hdiff;

	colormap;

	bp = &boxes[num];

	ldiff = bp->low[REDI] - bp->mean[REDI];
	ldiff *= ldiff;
	hdiff = bp->high[REDI] - bp->mean[REDI];
	hdiff *= hdiff;
	dist = MAX(ldiff, hdiff);

	ldiff = bp->low[GREENI] - bp->mean[GREENI];
	ldiff *= ldiff;
	hdiff = bp->high[GREENI] - bp->mean[GREENI];
	hdiff *= hdiff;
	dist += MAX(ldiff, hdiff);

	ldiff = bp->low[BLUEI] - bp->mean[BLUEI];
	ldiff *= ldiff;
	hdiff = bp->high[BLUEI] - bp->mean[BLUEI];
	hdiff *= hdiff;
	dist += MAX(ldiff, hdiff);

	dist = (float)sqrt((double)dist);
	/*
	 * Loop over all colors in the colormap, the ith entry of which
	 * corresponds to the ith box.
	 *
	 * If the centroid of a box is as close to any corner of the
	 * current box as is the centroid of the current box, add that
	 * box to the list of "neighbors" of the current box.
	 */
	HighR = (float)bp->high[REDI] + dist;
	HighG = (float)bp->high[GREENI] + dist;
	HighB = (float)bp->high[BLUEI] + dist;
	LowR = (float)bp->low[REDI] - dist;
	LowG = (float)bp->low[GREENI] - dist;
	LowB = (float)bp->low[BLUEI] - dist;
	for (i = j = 0, bp = boxes; i < colors; i++, bp++) {
		if (LowR <= bp->mean[REDI] && HighR >= bp->mean[REDI] &&
		    LowG <= bp->mean[GREENI] && HighG >= bp->mean[GREENI] &&
		    LowB <= bp->mean[BLUEI] && HighB >= bp->mean[BLUEI])
			neighbors[j++] = i;
	}

	return j;	/* Return the number of neighbors found. */
}

/*
 * Assign representative colors to every pixel in a given box through
 * the construction of the NearestColor array.  For each color in the
 * given box, we look at the list of neighbors passed to find the
 * one whose centroid is closest to the current color.
 */
static void makenearest(Box *boxes, int boxnum, int  nneighbors, int *neighbors, unsigned char *rgbmap, int bits, unsigned char *colormap[3],
int otherimages)
{
	register int n, b, g, r;
	double rdist, gdist, bdist, dist, mindist;
	int which, *np;
	Box *box;

	colormap;
	
	box = &boxes[boxnum];

	for (r = box->low[REDI]; r < box->high[REDI]; r++) {
		for (g = box->low[GREENI]; g < box->high[GREENI]; g++) {
			for (b = box->low[BLUEI]; b < box->high[BLUEI]; b++) {
/*
 * The following "if" is to avoid doing extra work when the RGBmap is
 * not going to be used for other images.
 */
			        if ((!otherimages) && 
				    (gHistogram[(((r<<bits)|g)<<bits)|b] == 0))
					continue;

				mindist = HUGE;
				/*
				 * Find the colormap entry which is
				 * closest to the current color.
				 */
				np = neighbors;
				for (n = 0; n < nneighbors; n++, np++) {
					rdist = r-boxes[*np].mean[REDI];
					gdist = g-boxes[*np].mean[GREENI];
					bdist = b-boxes[*np].mean[BLUEI];
					dist = rdist*rdist + gdist*gdist + bdist*bdist;
					if (dist < mindist) {
						mindist = dist;
						which = *np; 
					}
				}
				/*
				 * The colormap entry closest to this
				 * color is used as a representative.
				 */
				rgbmap[(((r<<bits)|g)<<bits)|b] = which;
			}
		}
	}
}
#endif /* slow_color_map */
