/* ppmquant.c - quantize the colors in a pixmap down to a specified number
**
** Copyright (C) 1989, 1991, 1993 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

#include "ppm.h"
#include "ppmcmap.h"

/* This define controls how big the histogram can get before the program
** rescales and tries again.
*/
/* #define MAXCOLORS 32767 */
#define MAXCOLORS 50000

/* This define controls how much the Floyd-Steinberg error diffusion
** can displace a color component.  The problem is that if you allow
** arbitrarily large displacements, the colors get too far away from
** the quantized colormap, and so a good matching pixel cannot be found
** and you get an ugly "singleton" pixel.  So instead we limit the error
** so that pixels can only be tweaked so far.  The number is relative
** to maxval, so for a maxval of 255, a limit of 0.125 means +- 31.
*/
#define FS_LIMIT 0.125

/* These defines control the algorithm used to decide which is the
** largest dimension of a box.
*/
/* #define LARGE_RGB */
#define LARGE_LUM

/* These defines control the algorithm used when choosing a color to
** represent a box.
*/
/* #define REP_CENTER_BOX */
/* #define REP_AVERAGE_COLORS */
#define REP_AVERAGE_PIXELS

/* These defines control how far the recombinant array code looks when
** it's trying to find a sharable array entry, which determines how much
** memory gets used.  The first define is used at initialization, for
** turning a user-supplied colormap into the recombinant array - it's not
** used at all if using median-cut to make the colormap.  The second
** define is used when doing Floyd-Steinberg error propagation, to add
** slightly-changed entries to the recombinant array.
*/
#define USERMAP_NEIGHBOR_DIST 2
#define FS_NEIGHBOR_DIST 5
/* Here's some data for a sample 200x450 image - for various values
** of the neighbor defines, the table shows how much memory gets used:
**     MB                         FS
**           0     1     2     3     4     5     6     7
**     U 0  ---  5.76  3.61  2.77  2.36  2.04  1.83  1.68
**     S 1  ---  4.40  2.65  2.01  1.59  1.39  1.24  1.15
**     E 2  ---  3.73  1.99  1.49  1.23  1.07  .960  .867
**     R 3  ---  3.26  1.77  1.29  1.10  .944  .841  .774
**     M 4  ---  3.24  1.73  1.22  .970  .849  .737  .714
**     A 5 10.5  2.97  1.45  1.14  .906  .787  .725  .633
**     P 6 10.7  3.01  1.54  1.17  .911  .833  .690  .614
** The other constraints besides memory use are collisions and search time.
** If USERMAP_NEIGHBOR_DIST is too large, shared array entries collide
** and have to be split, and for each collision there's a small chance that
** a colormap entry will get lost.  We want to avoid that, so we set
** USERMAP_NEIGHBOR_DIST to two.  For the FS_NEIGHBOR_DIST, it looks like
** the larger it gets the less memory is used; however, the larger numbers
** mean more time spent searching for array entries that can be shared.
** Plus, unsplittable collisions can happen here too.  Four or five looks
** like a good compromise.
*/

/* If maxval is 127, a complete 3-D unshared array is only two megabytes,
** and at maxval 63 it's only a quarter meg.  This define sets the
** threshold below which we don't bother with sharing arrays.
**
** Note that at some point in the not too distant future, it will become
** reasonable to allocate a full 256*256*256 array - that's only 16 MB.
** At that point, the recombinant array code will be obsolete.
*/
#define COMPLETE_MAXVAL 99

typedef enum { NotSorted, RedSorted, GreenSorted, BlueSorted } sort_type;

typedef struct box* box_vector;
struct box
    {
    int ind;
    int colors;
    long sum;
    pixel color;
    sort_type sorted_by;
    pixval minr, maxr, ming, maxg, minb, maxb;
    };

static box_vector median_cut( colorhist_vector chv, int colors, long sum, pixval maxval, int newcolors );
static void choose_reps( box_vector bv, int newcolors, colorhist_vector chv, pixval maxval );
static void bound_box( struct box* boxP, colorhist_vector chv );
static void add_box_to_array( pixel**** color_array, struct box* box, pixval maxval );
static void add_pixel_to_array( pixel**** color_array, int neighbor_dist, int r, int g, int b, pixel* newpP, pixval maxval );
static void map_colors( pixel **pixels, pixel* pixelrow, int readlines, FILE* ifp, int rows, int cols, pixval maxval, int format, pixel**** color_array, colorhist_vector chv, int newcolors, int floyd, int fs_limit, int fs_neighbor_dist );
static int find_best_color( pixel* pP, pixval maxval, colorhist_vector chv, int newcolors );


int
main( int argc, char* argv[] )
    {
    FILE* ifp;
    pixel** pixels = (pixel**) 0;
    pixel** mappixels;
    pixel* pixelrow = (pixel*) 0;
    int readlines;
    int argn, rows, cols, format = -1, maprows, mapcols, i;
    pixval mapmaxval, maxval;
    int newcolors, colors;
    colorhist_vector chv;
    int floyd, fs_limit, writemap;
    int usermap_neighbor_dist, fs_neighbor_dist;
    pixel**** color_array;
    char* usage = "[-nofloyd|-nofs] [-nolimit] <ncolors> [ppmfile]\n                 [-nofloyd|-nofs] [-nolimit] -map mapfile [ppmfile]";

    ppm_init( &argc, argv );

    argn = 1;
    floyd = 1;
    fs_limit = 1;
    writemap = 0;
    mappixels = (pixel**) 0;

    while ( argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0' )
	{
	if ( pm_keymatch( argv[argn], "-fs", 2 ) ||
	     pm_keymatch( argv[argn], "-floyd", 2 ) )
	    floyd = 1;
	else if ( pm_keymatch( argv[argn], "-nofs", 4 ) ||
		  pm_keymatch( argv[argn], "-nofloyd", 4 ) )
	    floyd = 0;
	else if ( pm_keymatch( argv[argn], "-limit", 2 ) )
	    fs_limit = 1;
	else if ( pm_keymatch( argv[argn], "-nolimit", 4 ) )
	    fs_limit = 0;
	else if ( pm_keymatch( argv[argn], "-writemap", 2 ) )
	    writemap = 1;
	else if ( pm_keymatch( argv[argn], "-map", 2 ) )
	    {
	    ++argn;
	    if ( argn == argc )
		pm_usage( usage );
	    ifp = pm_openr( argv[argn] );
	    mappixels = ppm_readppm( ifp, &mapcols, &maprows, &mapmaxval );
	    pm_closer( ifp );
	    if ( mapcols == 0 || maprows == 0 )
		pm_error( "null colormap??" );
	    }
	else
	    pm_usage( usage );
	++argn;
	}

    if ( mappixels == (pixel**) 0 )
	{
	if ( argn == argc )
	    pm_usage( usage );
	if ( sscanf( argv[argn], "%d", &newcolors ) != 1 )
	    pm_usage( usage );
	if ( newcolors < 1 )
	    pm_error( "number of colors must be >= 1" );
	++argn;
	}

    if ( argn != argc )
	{
	ifp = pm_openr( argv[argn] );
	++argn;
	}
    else
	ifp = stdin;

    if ( argn != argc )
	pm_usage( usage );

    if ( mappixels == (pixel**) 0 )
	{
	box_vector bv;

	/* Read in the image. */
	pixels = ppm_readppm( ifp, &cols, &rows, &maxval );
	pm_closer( ifp );
	readlines = 0;

	/* Attempt to make a histogram of the colors, unclustered.
	** If at first we don't succeed, lower maxval to increase color
	** coherence and try again.  This will eventually terminate, with
	** maxval at worst 15, since 32^3 is approximately MAXCOLORS.
	*/
	for ( ; ; )
	    {
	    register int rmv, newmaxval;
	    register int row, col;
	    register pixel* pP;

	    pm_message( "making histogram..." );
	    chv = ppm_computecolorhist(
		pixels, cols, rows, MAXCOLORS, &colors );
	    if ( chv != (colorhist_vector) 0 )
		break;
	    pm_message( "too many colors!" );
	    rmv = maxval;
	    newmaxval = rmv / 2;
	    pm_message(
	  "scaling colors from maxval=%d to maxval=%d to improve clustering...",
		    rmv, newmaxval );
	    for ( row = 0; row < rows; ++row )
		for ( col = 0, pP = pixels[row]; col < cols; ++col, ++pP )
		    PPM_DEPTH( *pP, *pP, rmv, newmaxval );
	    maxval = newmaxval;
	    }
	pm_message(
	    "%d different color%s found", colors, (colors == 1) ? "" : "s" );

	/* Possible short circuit. */
	if ( colors <= newcolors )
	    {
	    pm_message( "few enough colors now; copying image..." );
	    ppm_writeppm( stdout, pixels, cols, rows, maxval, 0 );
	    pm_closew( stdout );
	    exit( 0 );
	    }

	/* Apply median-cut to the histogram, making the new colormap. */
	pm_message(
	    "choosing %d color%s...", newcolors, (newcolors == 1) ? "" : "s" );
	bv = median_cut( chv, colors, rows * cols, maxval, newcolors );

	/* Choose representative colors for the boxes. */
	choose_reps( bv, newcolors, chv, maxval );

	ppm_freecolorhist( chv );

	/* Short-circuit is we're writing just the map. */
	if ( writemap )
	    {
	    pixel** mp;

	    mp = ppm_allocarray( newcolors, 1 );
	    for ( i = 0; i < newcolors; ++i )
		mp[0][i] = bv[i].color;
	    ppm_writeppm( stdout, mp, newcolors, 1, maxval, 0 );
	    exit( 0 );
	    }

	/* Initialize recombinant array for color mapping. */
	color_array = (pixel****) pm_allocrow(
	    (int) maxval + 1, sizeof(pixel***) );
	for ( i = 0; i <= maxval; ++i )
	    color_array[i] = (pixel***) 0;
	
	/* Convert the box list returned by median_cut() into a 3-D
	** recombinant array of pixel*, indexed by [r][g][b].
	*/
	for ( i = 0; i < newcolors; ++i )
	    add_box_to_array( color_array, &(bv[i]), maxval );

	/* And also convert the box list into a color list, for use
	** by find_best_color().
	*/
	chv = (colorhist_vector) pm_allocrow(
	    newcolors, sizeof(struct colorhist_item) );
	for ( i = 0; i < newcolors; ++i )
	    chv[i].color = bv[i].color;

	free( (char*) bv );
	}
    else
	{
	register pixel* pP;

	/* Alternate path for user-supplied colormap.  Read in the
	** image header.
	*/
	ppm_readppminit( ifp, &cols, &rows, &maxval, &format );
	pixelrow = ppm_allocrow( cols );
	readlines = 1;

	/* Turn mappixels into a colormap. */
	if ( mapmaxval != maxval )
	    {
	    register int row, col;

	    if ( mapmaxval > maxval )
		pm_message( "rescaling colormap colors" );
	    for ( row = 0; row < maprows; ++row )
		for ( col = 0, pP = mappixels[row]; col < mapcols; ++col, ++pP )
		    PPM_DEPTH( *pP, *pP, mapmaxval, maxval );
	    }
	chv = ppm_computecolorhist(
	    mappixels, mapcols, maprows, MAXCOLORS, &newcolors );
	if ( chv == (colorhist_vector) 0 )
	    pm_error( "too many colors in colormap!" );
	ppm_freearray( mappixels, maprows );
	pm_message(
	    "%d different color%s found in colormap", newcolors,
	    (newcolors == 1) ? "" : "s" );

	/* Initialize recombinant array for color mapping. */
	color_array = (pixel****) pm_allocrow(
	    (int) maxval + 1, sizeof(pixel***) );
	for ( i = 0; i <= maxval; ++i )
	    color_array[i] = (pixel***) 0;

	if ( maxval <= COMPLETE_MAXVAL )
	    usermap_neighbor_dist = 0;
	else
	    usermap_neighbor_dist = USERMAP_NEIGHBOR_DIST;

	/* Convert the colormap into a 3-D recombinant array of pixel*,
	** indexed by [r][g][b].
	*/
	for ( i = 0; i < newcolors; ++i )
	    {
	    pP = &(chv[i].color);
	    add_pixel_to_array(
		color_array, usermap_neighbor_dist, 
		PPM_GETR( *pP ), PPM_GETG( *pP ), PPM_GETB( *pP ),
		pP, maxval );
	    }
	}

    if ( maxval <= COMPLETE_MAXVAL )
	fs_neighbor_dist = 0;
    else
	fs_neighbor_dist = FS_NEIGHBOR_DIST;
    
    /* Map to the new colors and write the file. */
    pm_message( "mapping image to new colors..." );
    map_colors(
	pixels, pixelrow, readlines, ifp, rows, cols, maxval, format,
	color_array, chv, newcolors, floyd, fs_limit, fs_neighbor_dist );

    /* All done. */
    pm_closew( stdout );
    exit( 0 );
    }


/* qsort comparison functions */

static int
red_compare( const void* v1, const void* v2 )
    {
    const colorhist_vector ch1 = (const colorhist_vector) v1;
    const colorhist_vector ch2 = (const colorhist_vector) v2;

    return (int) PPM_GETR( ch1->color ) - (int) PPM_GETR( ch2->color );
    }

static int
green_compare( const void* v1, const void* v2 )
    {
    const colorhist_vector ch1 = (const colorhist_vector) v1;
    const colorhist_vector ch2 = (const colorhist_vector) v2;

    return (int) PPM_GETG( ch1->color ) - (int) PPM_GETG( ch2->color );
    }

static int
blue_compare( const void* v1, const void* v2 )
    {
    const colorhist_vector ch1 = (const colorhist_vector) v1;
    const colorhist_vector ch2 = (const colorhist_vector) v2;

    return (int) PPM_GETB( ch1->color ) - (int) PPM_GETB( ch2->color );
    }


/* Here is the fun part, the median-cut colormap generator.  This is based
** on Paul Heckbert's paper "Color Image Quantization for Frame Buffer
** Display", SIGGRAPH '82 Proceedings, page 297.
*/
static box_vector
median_cut( colorhist_vector chv, int colors, long sum, pixval maxval, int newcolors )
    {
    box_vector bv;
    register int bi, i, j;
    int pi;
    int boxes;
    struct box newbox, tempbox;
    int num_to_move;

    bv = (box_vector) pm_allocrow( newcolors, sizeof(struct box) );

    /* Set up the initial box. */
    bv[0].ind = 0;
    bv[0].colors = colors;
    bv[0].sum = sum;
    bv[0].sorted_by = NotSorted;
    bound_box( &bv[0], chv );
    boxes = 1;

    /* Main median-cut loop: split boxes until we have enough. */
    while ( boxes < newcolors )
	{
	register int indx, clrs;
	long sm, halfsum, lowersum, plowersum;

	/* Find the first splittable box. */
	for ( bi = 0; bi < boxes; ++bi )
	    if ( bv[bi].colors >= 2 )
		break;
	if ( bi == boxes )
	    break;	/* ran out of colors! */
	indx = bv[bi].ind;
	clrs = bv[bi].colors;
	sm = bv[bi].sum;

	/* Find the largest splittable dimension, and sort by that
	** component.  I have included two methods for determining the
	** largest dimension: first by simply comparing the range in RGB
	** space, and second by transforming into luminosities before the
	** comparison.  You can switch which method is used by switching
	** the commenting on the LARGE_ defines at the beginning of this
	** source file.
	*/
#ifdef LARGE_RGB
	{
	int rd, gd, bd;

	rd = (int) bv[bi].maxr - (int) bv[bi].minr;
	gd = (int) bv[bi].maxg - (int) bv[bi].ming;
	bd = (int) bv[bi].maxb - (int) bv[bi].minb;

	if ( rd >= gd && rd >= bd && rd > 0 )
	    {
	    if ( bv[bi].sorted_by != RedSorted )
		{
		qsort(
		    (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
		    red_compare );
		bv[bi].sorted_by = RedSorted;
		}
	    }
	else if ( gd >= bd && gd > 0 )
	    {
	    if ( bv[bi].sorted_by != GreenSorted )
		{
		qsort(
		    (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
		    green_compare );
		bv[bi].sorted_by = GreenSorted;
		}
	    }
	else if ( bd > 0 )
	    {
	    if ( bv[bi].sorted_by != BlueSorted )
		{
		qsort(
		    (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
		    blue_compare );
		bv[bi].sorted_by = BlueSorted;
		}
	    }
	else
	    pm_error( "largest box isn't splittable?!?" );
	}
#endif /*LARGE_RGB*/
#ifdef LARGE_LUM
	{
	pixel p;
	int rd, gd, bd;
	float rl, gl, bl;

	rd = (int) bv[bi].maxr - (int) bv[bi].minr;
	gd = (int) bv[bi].maxg - (int) bv[bi].ming;
	bd = (int) bv[bi].maxb - (int) bv[bi].minb;
	PPM_ASSIGN( p, rd, 0, 0 );
	rl = PPM_LUMIN(p);
	PPM_ASSIGN( p, 0, gd, 0 );
	gl = PPM_LUMIN(p);
	PPM_ASSIGN( p, 0, 0, bd );
	bl = PPM_LUMIN(p);

	if ( rl >= gl && rl >= bl && rd > 0 )
	    {
	    if ( bv[bi].sorted_by != RedSorted )
		{
		qsort(
		    (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
		    red_compare );
		bv[bi].sorted_by = RedSorted;
		}
	    }
	else if ( gl >= bl && gd > 0 )
	    {
	    if ( bv[bi].sorted_by != GreenSorted )
		{
		qsort(
		    (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
		    green_compare );
		bv[bi].sorted_by = GreenSorted;
		}
	    }
	else if ( bd > 0 )
	    {
	    if ( bv[bi].sorted_by != BlueSorted )
		{
		qsort(
		    (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
		    blue_compare );
		bv[bi].sorted_by = BlueSorted;
		}
	    }
	else
	    pm_error( "largest box isn't splittable?!?" );
	}
#endif /*LARGE_LUM*/
	
	/* Find the legal split point nearest the median based on counts,
	** so that about half the pixels (not colors, pixels) are in each
	** part.  A legal split point is one where the color component we
	** just sorted on changes value.  We have three copies of the
	** loop to avoid procedure calls / branches inside the loop.
	*/
	halfsum = sm / 2;
	lowersum = chv[indx].value;
	switch ( bv[bi].sorted_by )
	    {
	    case RedSorted:
	    for ( i = 1, pi = -1; i < clrs; ++i )
		{
		if ( PPM_GETR( chv[indx+i-1].color ) !=
		     PPM_GETR( chv[indx+i].color ) )
		    {
		    if ( lowersum >= halfsum )
			break;
		    pi = i;
		    plowersum = lowersum;
		    }
		lowersum += chv[indx + i].value;
		}
	    break;

	    case GreenSorted:
	    for ( i = 1, pi = -1; i < clrs; ++i )
		{
		if ( PPM_GETG( chv[indx+i-1].color ) !=
		     PPM_GETG( chv[indx+i].color ) )
		    {
		    if ( lowersum >= halfsum )
			break;
		    pi = i;
		    plowersum = lowersum;
		    }
		lowersum += chv[indx + i].value;
		}
	    break;

	    case BlueSorted:
	    for ( i = 1, pi = -1; i < clrs; ++i )
		{
		if ( PPM_GETB( chv[indx+i-1].color ) !=
		     PPM_GETB( chv[indx+i].color ) )
		    {
		    if ( lowersum >= halfsum )
			break;
		    pi = i;
		    plowersum = lowersum;
		    }
		lowersum += chv[indx + i].value;
		}
	    break;

	    default:
	    pm_error( "unknown sort type?!?" );
	    }

	if ( pi == -1 )
	    { /* No previous split point, use current. */
	    if ( i == clrs )
		/* No previous or current split point? */
		pm_error( "couldn't find a split point?!?" );
	    }
	else if ( i == clrs )
	    { /* No current split point, use previous. */
	    i = pi;
	    lowersum = plowersum;
	    }
	else
	    { /* Have both previous and current split points, use best. */
	    if ( halfsum - plowersum < lowersum - halfsum )
		{ /* Previous is better. */
		i = pi;
		lowersum = plowersum;
		}
	    }

	/* Split the box between i-1 and i. */
	bv[bi].colors = i;
	bv[bi].sum = lowersum;
	/* Recompute bounding box instead of just splitting on the sorted
	** component, so we can shrink boxes better.
	*/
	bound_box( &bv[bi], chv );

	newbox.ind = indx + i;
	newbox.colors = clrs - i;
	newbox.sum = sm - lowersum;
	newbox.sorted_by = bv[bi].sorted_by;
	/* And recompute instead of splitting here too. */
	bound_box( &newbox, chv );

	/* Sort to bring the biggest box to the top.  We don't need to
	** resort the whole box list, only the first box and the new box
	** are out of order.
	*/
	for ( i = bi + 1; i < boxes && lowersum < bv[i].sum; ++i)
	    ;
	num_to_move = ( i - 1 ) - bi;
	if ( num_to_move > 0 )
	    {
	    tempbox = bv[bi];
	    for ( j = 0; j < num_to_move; ++j )
		bv[bi + j] = bv[bi + j + 1];
	    bv[i - 1] = tempbox;
	    }

	/* Now, insert the new box */
	for ( i = bi; i < boxes && newbox.sum <= bv[i].sum; ++i)
	    ;
	num_to_move = boxes - i;
	if ( num_to_move > 0 )
	    for ( j = num_to_move - 1; j >= 0; --j )
		bv[i + j + 1] = bv[i + j];
	bv[i] = newbox;

	++boxes;
	}

    /* Done selecting colors. */
    return bv;
    }


/* Choose a representative color for each box.  There are a number of
** possible ways to make this choice.  One would be to choose the center
** of the box; this ignores any structure within the boxes.  Another
** method would be to average all the colors in the box - this is the
** method specified in Heckbert's paper.  A third method is to average
** all the pixels in the box (i.e. a color-average weighted by the
** prevalence of each color).  You can switch which method is used by
** switching the commenting on the REP_ defines at the beginning of this
** source file.
*/
static void
choose_reps( box_vector bv, int newcolors, colorhist_vector chv, pixval maxval )
    {
    register int ind, colors;
    register long r, g, b, sum;
    register int i, j;

    for ( i = 0; i < newcolors; ++i )
	{
#ifdef REP_CENTER_BOX
	PPM_ASSIGN(
	    bv[i].color, ( (int) bv[i].minr + (int) bv[i].maxr ) / 2,
	    ( (int) bv[i].ming + (int) bv[i].maxg ) / 2,
	    ( (int) bv[i].minb + (int) bv[i].maxb ) / 2 );
#endif /*REP_CENTER_BOX*/
#ifdef REP_AVERAGE_COLORS
	ind = bv[i].ind;
	colors = bv[i].colors;
	r = g = b = 0;

	for ( j = 0; j < colors; ++j )
	    {
	    r += PPM_GETR( chv[ind + j].color );
	    g += PPM_GETG( chv[ind + j].color );
	    b += PPM_GETB( chv[ind + j].color );
	    }
	r = r / colors;
	g = g / colors;
	b = b / colors;
	PPM_ASSIGN( bv[i].color, r, g, b );
#endif /*REP_AVERAGE_COLORS*/
#ifdef REP_AVERAGE_PIXELS
	ind = bv[i].ind;
	colors = bv[i].colors;
	sum = bv[i].sum;
	r = g = b = 0;

	for ( j = 0; j < colors; ++j )
	    {
	    r += PPM_GETR( chv[ind + j].color ) * chv[ind + j].value;
	    g += PPM_GETG( chv[ind + j].color ) * chv[ind + j].value;
	    b += PPM_GETB( chv[ind + j].color ) * chv[ind + j].value;
	    }
	r = r / sum;
	if ( r > maxval ) r = maxval;	/* avoid math errors */
	g = g / sum;
	if ( g > maxval ) g = maxval;
	b = b / sum;
	if ( b > maxval ) b = maxval;
	PPM_ASSIGN( bv[i].color, r, g, b );
#endif /*REP_AVERAGE_PIXELS*/
	}
    }


/* Find the minimum and maximum of each color component in a box - the
** boundaries.
*/
static void
bound_box( boxP, chv )
    struct box* boxP;
    colorhist_vector chv;
    {
    register int ind, colors, i;
    register pixval minr, maxr, ming, maxg, minb, maxb, v;

    ind = boxP->ind;
    colors = boxP->colors;
    minr = maxr = PPM_GETR( chv[ind].color );
    ming = maxg = PPM_GETG( chv[ind].color );
    minb = maxb = PPM_GETB( chv[ind].color );
    for ( i = 1; i < colors; ++i )
	{
	v = PPM_GETR( chv[ind + i].color );
	if ( v < minr ) minr = v;
	else if ( v > maxr ) maxr = v;
	v = PPM_GETG( chv[ind + i].color );
	if ( v < ming ) ming = v;
	else if ( v > maxg ) maxg = v;
	v = PPM_GETB( chv[ind + i].color );
	if ( v < minb ) minb = v;
	else if ( v > maxb ) maxb = v;
	}
    boxP->minr = minr;
    boxP->maxr = maxr;
    boxP->ming = ming;
    boxP->maxg = maxg;
    boxP->minb = minb;
    boxP->maxb = maxb;
    }


/* Convert a box-vector entry into a 3-D recombinant array of pixel*,
** indexed by [r][g][b].
**
** This is used to convert a set of boxes produced by median-cut into
** the initial recombinant array.
*/
static void
add_box_to_array( pixel**** color_array, struct box* box, pixval maxval )
    {
    int r, g, b;
    pixel*** blankredvec;
    pixel** blankgreenvec;
    int blankred_used, blankgreen_used, newp_used;
    pixel*** redvec;
    pixel** greenvec;
    pixel* newpP;
    pixel* pP;
    pixel*** newredvec;
    pixel** newgreenvec;
    pixel*** prevredvec;
    pixel** prevgreenvec;

    /* Create two blank intermediate vectors and the bottom-level pixel. */
    blankredvec = (pixel***) pm_allocrow( (int) maxval + 1, sizeof(pixel**) );
    for ( g = 0; g <= maxval; ++g )
	blankredvec[g] = (pixel**) 0;
    blankred_used = 0;

    blankgreenvec = (pixel**) pm_allocrow( (int) maxval + 1, sizeof(pixel*) );
    for ( b = 0; b <= maxval; ++b )
	blankgreenvec[b] = (pixel*) 0;
    blankgreen_used = 0;

    newpP = (pixel*) malloc( sizeof(pixel) );
    if ( newpP == (pixel*) 0 )
	pm_error( "out of memory" );
    *newpP = box->color;
    newp_used = 0;

    /* Fill in any blank red entries, copy any existing ones. */
    for ( prevredvec = (pixel***) 0, r = box->minr; r <= box->maxr; ++r )
	{
	redvec = color_array[r];
	if ( redvec == (pixel***) 0 )
	    {
	    color_array[r] = blankredvec;
	    blankred_used = 1;
	    }
	else
	    {
	    if ( redvec != prevredvec )
		{
		newredvec = (pixel***) pm_allocrow(
		    (int) maxval + 1, sizeof(pixel**) );
		for ( g = 0; g <= maxval; ++g )
		    newredvec[g] = redvec[g];
		prevredvec = redvec;
		}
	    color_array[r] = newredvec;
	    }
	}

    /* Now loop through unique reds and handle next level down. */
    for ( prevredvec = (pixel***) 0, r = box->minr; r <= box->maxr; ++r )
	{
	redvec = color_array[r];
	if ( redvec != prevredvec )
	    {
	    prevredvec = redvec;

	    /* Fill in any blank green entries, copy any existing ones. */
	    for ( prevgreenvec = (pixel**) 0, g = box->ming;
		  g <= box->maxg; ++g )
		{
		greenvec = redvec[g];
		if ( greenvec == (pixel**) 0 )
		    {
		    redvec[g] = blankgreenvec;
		    blankgreen_used = 1;
		    }
		else
		    {
		    if ( greenvec != prevgreenvec )
			{
			newgreenvec = (pixel**) pm_allocrow(
			    (int) maxval + 1, sizeof(pixel*) );
			for ( b = 0; b <= maxval; ++b )
			    newgreenvec[b] = greenvec[b];
			prevgreenvec = greenvec;
			}
		    redvec[g] = newgreenvec;
		    }
		}

	    /* Now loop through unique greens and handle next level down. */
	    for ( prevgreenvec = (pixel**) 0, g = box->ming;
		  g <= box->maxg; ++g )
		{
		greenvec = redvec[g];
		if ( greenvec != prevgreenvec )
		    {
		    prevgreenvec = greenvec;

		    /* Fill in any missing blue entries. */
		    for ( b = box->minb; b <= box->maxb; ++b )
			{
			pP = greenvec[b];
			if ( pP == (pixel*) 0 )
			    {
			    greenvec[b] = newpP;
			    newp_used = 1;
			    }
			else if ( pP == newpP )
			    break;  /* we've already done this one */
			else if ( ! PPM_EQUAL( *pP, *newpP ) )
			    /* Oops.  This should never happen, because
			    ** the median-cut box splitter should never
			    ** produce overlapping boxes.
			    */
			    pm_error( "overlapping non-equal color boxes?!?" );
			}
		    }
		}
	    }
	}

    if ( ! blankred_used )
	free( (char*) blankredvec );
    if ( ! blankgreen_used )
	free( (char*) blankgreenvec );
    if ( ! newp_used )
	free( (char*) newpP );
    }


/* Convert a colormap entry into a 3-D recombinant array of pixel*, indexed
** by [r][g][b].
**
** This is used both to convert a user-supplied colormap into the initial
** recombinant array, and during Floyd-Steinberg error propagation to add
** slightly modified entries to the array.
*/
static void
add_pixel_to_array( pixel**** color_array, int neighbor_dist, int r, int g, int b, pixel* newpP, pixval maxval )
    {
    register int i;
    register pixel*** redvec;
    register pixel** greenvec;
    register pixel* pP;

    /* Fill in the red entry, if missing. */
    redvec = color_array[r];
    if ( redvec == (pixel***) 0 )
	{
	/* Check red neighbors to see if we can use them. */
	for ( i = 1; i <= neighbor_dist; ++i )
	    {
	    if ( r-i >= 0 && color_array[r-i] != (pixel***) 0 &&
		 ( color_array[r-i][g] == (pixel**) 0 ||
		   color_array[r-i][g][b] == (pixel*) 0 ||
		   PPM_EQUAL( *color_array[r-i][g][b], *newpP ) ) )
		{
		redvec = color_array[r] = color_array[r-i];
		break;
		}
	    else if ( r+i <= (int) maxval && color_array[r+i] != (pixel***) 0 &&
		      ( color_array[r+i][g] == (pixel**) 0 ||
			color_array[r+i][g][b] == (pixel*) 0 ||
			PPM_EQUAL( *color_array[r+i][g][b], *newpP ) ) )
		{
		redvec = color_array[r] = color_array[r+i];
		break;
		}
	    }
	if ( redvec == (pixel***) 0 )
	    {
	    /* Nope, gotta make a new red entry. */
	    redvec = color_array[r] = (pixel***) pm_allocrow(
		(int) maxval + 1, sizeof(pixel**) );
	    for ( i = 0; i <= maxval; ++i )
		redvec[i] = (pixel**) 0;
	    }
	}

    /* Fill in the green entry, if missing. */
    greenvec = redvec[g];
    if ( greenvec == (pixel**) 0 )
	{
	/* Check green neighbors to see if we can use them. */
	for ( i = 1; i <= neighbor_dist; ++i )
	    {
	    if ( g-i >= 0 && redvec[g-i] != (pixel**) 0 &&
		 ( redvec[g-i][b] == (pixel*) 0 ||
		   PPM_EQUAL( *redvec[g-i][b], *newpP ) ) )
		{
		greenvec = redvec[g] = redvec[g-i];
		break;
		}
	    else if ( g+i <= (int) maxval && redvec[g+i] != (pixel**) 0 &&
		      ( redvec[g+i][b] == (pixel*) 0 ||
			PPM_EQUAL( *redvec[g+i][b], *newpP ) ) )
		{
		greenvec = redvec[g] = redvec[g+i];
		break;
		}
	    }
	if ( greenvec == (pixel**) 0 )
	    {
	    /* Nope, gotta make a new green entry. */
	    greenvec = redvec[g] = (pixel**) pm_allocrow(
		(int) maxval + 1, sizeof(pixel*) );
	    for ( i = 0; i <= maxval; ++i )
		greenvec[i] = (pixel*) 0;
	    }
	}

    /* Fill in the blue entry. */
    pP = greenvec[b];
    if ( pP == (pixel*) 0 )
	{
	/* Check blue neighbors to see if we can use them. */
	for ( i = 1; i <= neighbor_dist; ++i )
	    {
	    if ( b-i >= 0 && greenvec[b-i] != (pixel*) 0 &&
		 PPM_EQUAL( *greenvec[b-i], *newpP ) )
		{
		greenvec[b] = greenvec[b-i];
		break;
		}
	    else if ( b+i <= (int) maxval && greenvec[b+i] != (pixel*) 0 &&
		      PPM_EQUAL( *greenvec[b+i], *newpP ) )
		{
		greenvec[b] = greenvec[b+i];
		break;
		}
	    }
	if ( pP == (pixel*) 0 )
	    {
	    /* Nope, gotta make a new blue entry. */
	    pP = greenvec[b] = (pixel*) malloc( sizeof(pixel) );
	    if ( pP == (pixel*) 0 )
		pm_error( "out of memory" );
	    *pP = *newpP;
	    }
	}
    else
	{
	if ( ! PPM_EQUAL( *pP, *newpP ) )
	    {
	    /* Oops, a collision.  We have to split off a clean set
	    ** of entries from the ones we oh-so-cleverly re-used before.
	    ** Them's the breaks, that's why this data structure is
	    ** so much fun.
	    **
	    ** Now, we only want to split off entries that are in fact
	    ** currently being shared - if they're not being shared there's
	    ** no point in making a new copy.  First check red entry.
	    */
	    for ( i = 0; i <= maxval; ++i )
		{
		if ( i != r && color_array[i] == color_array[r] )
		    {
		    /* Yes, it is - split it off. */
		    redvec = (pixel***) pm_allocrow(
			(int) maxval + 1, sizeof(pixel**) );
		    for ( i = 0; i <= maxval; ++i )
			redvec[i] = color_array[r][i];
		    color_array[r] = redvec;
		    }
		}
	    /* Now check whether the green entry is currently being shared. */
	    for ( i = 0; i <= maxval; ++i )
		{
		if ( i != g && color_array[r][i] == color_array[r][g] )
		    {
		    /* Yes, it is - split it off. */
		    greenvec = (pixel**) pm_allocrow(
			(int) maxval + 1, sizeof(pixel*) );
		    for ( i = 0; i <= maxval; ++i )
			greenvec[i] = color_array[r][g][i];
		    color_array[r][g] = greenvec;
		    }
		}
	    /* And check whether the blue entry is currently being shared. */
	    for ( i = 0; i <= maxval; ++i )
		{
		if ( i != b && color_array[r][g][i] == color_array[r][g][b] )
		    {
		    /* Yes, it is - split it off. */
		    pP = (pixel*) malloc( sizeof(pixel) );
		    if ( pP == (pixel*) 0 )
			pm_error( "out of memory" );
		    color_array[r][g][b] = pP;
		    }
		}
	    /* Finally, we can set the new pixel value.
	    **
	    ** It's possible for none of the above three split attempts
	    ** to succeed.  That can happen if entries were already split
	    ** off before.  In that case, setting the new pixel value
	    ** might possibly be overwriting a real existing pixel value.
	    ** This is not a big deal for three reasons.  One, it happens
	    ** very rarely.  Two, when it does happen, most of the time
	    ** the old value is not real, it's just left over from when
	    ** the entry was shared.  And three, if it is in fact a real
	    ** value, the new value will be close to the old value, due
	    ** to the nature of the color quantization problem.
	    */
	    *color_array[r][g][b] = *newpP;
	    }
	}
    }


/* Write out the file using only the new colormap. */
#define FS_SCALE 1024
static void map_colors( pixel **pixels, pixel* pixelrow, int readlines, FILE* ifp, int rows, int cols, pixval maxval, int format, pixel**** color_array, colorhist_vector chv, int newcolors, int floyd, int fs_limit, int fs_neighbor_dist )
    {
    long *thisrerr;
    long *nextrerr;
    long *thisgerr;
    long *nextgerr;
    long *thisberr;
    long *nextberr;
    int row, col;
    int fs_direction;
    int limitcol;
    pixval** limited_add_tab;
    register pixel* pP;
    long r, g, b, err;
    pixel*** redvec;
    pixel** greenvec;
    pixel* newpP;
    int i;

    ppm_writeppminit( stdout, cols, rows, maxval, 0 );

    if ( floyd )
	{
	/* Initialize Floyd-Steinberg error vectors. */
	thisrerr = (long*) pm_allocrow( cols + 2, sizeof(long) );
	nextrerr = (long*) pm_allocrow( cols + 2, sizeof(long) );
	thisgerr = (long*) pm_allocrow( cols + 2, sizeof(long) );
	nextgerr = (long*) pm_allocrow( cols + 2, sizeof(long) );
	thisberr = (long*) pm_allocrow( cols + 2, sizeof(long) );
	nextberr = (long*) pm_allocrow( cols + 2, sizeof(long) );
	pm_srandom( 0 );
	for ( col = 0; col < cols + 2; ++col )
	    {
	    thisrerr[col] = pm_random() % ( FS_SCALE * 2 ) - FS_SCALE;
	    thisgerr[col] = pm_random() % ( FS_SCALE * 2 ) - FS_SCALE;
	    thisberr[col] = pm_random() % ( FS_SCALE * 2 ) - FS_SCALE;
	    /* (random errors in [-1 .. 1], scaled) */
	    }
	fs_direction = 1;

	/* Initialize limited addition table.  For maxval=255, takes up 128K. */
	limited_add_tab = (pixval**) pm_allocarray(
	    maxval * 2 + 1, maxval + 1, sizeof(pixval) );
	b = maxval * FS_LIMIT;
	for ( r = 0; r <= maxval; ++r )
	    {
	    limited_add_tab[r] += maxval;
	    for ( err = - (long) maxval; err <= (long) maxval; ++err )
		{
		g = r + err;
		if ( g < 0 ) g = 0;
		else if ( g > maxval ) g = maxval;
		if ( fs_limit )
		    {
		    if ( g < r - b ) g = r - b;
		    else if ( g > r + b ) g = r + b;
		    }
		limited_add_tab[r][err] = g;
		}
	    }
	}

    for ( row = 0; row < rows; ++row )
	{
	if ( floyd )
	    for ( col = 0; col < cols + 2; ++col )
		nextrerr[col] = nextgerr[col] = nextberr[col] = 0;
	if ( readlines )
	    ppm_readppmrow( ifp, pixelrow, cols, maxval, format );
	else
	    pixelrow = pixels[row];
	if ( ( ! floyd ) || fs_direction )
	    {
	    col = 0;
	    limitcol = cols;
	    pP = pixelrow;
	    }
	else
	    {
	    col = cols - 1;
	    limitcol = -1;
	    pP = &(pixelrow[col]);
	    }

	do
	    {
	    r = PPM_GETR( *pP );
	    g = PPM_GETG( *pP );
	    b = PPM_GETB( *pP );
	    if ( floyd )
		{
		/* Use Floyd-Steinberg errors to adjust actual color. */
		r = limited_add_tab[r][thisrerr[col + 1] / FS_SCALE];
		g = limited_add_tab[g][thisgerr[col + 1] / FS_SCALE];
		b = limited_add_tab[b][thisberr[col + 1] / FS_SCALE];
		PPM_ASSIGN( *pP, r, g, b );
		}

	    /* Check closest-color recombinant array. */
	    if ( ( redvec = color_array[r] ) != (pixel***) 0 &&
		 ( greenvec = redvec[g] ) != (pixel**) 0 &&
		 ( newpP = greenvec[b] ) != (pixel*) 0 )
		{
		/* Change *pP to the closest color in the new colormap. */
		*pP = *newpP;
		}
	    else
		{
		/* Got to add a new entry. */
		i = find_best_color( pP, maxval, chv, newcolors );
		add_pixel_to_array(
		    color_array, fs_neighbor_dist, r, g, b,
		    &(chv[i].color), maxval );

		/* And change *pP to new closest color. */
		*pP = chv[i].color;
		}

	    if ( floyd )
		{
		/* Propagate Floyd-Steinberg error terms.  We optimize by
		** multiplying each error by FS_SCALE/16 instead of dividing
		** each sub-error by 16.  Also do decay factor here.
		*/
		if ( fs_direction )
		    {
		    err = ( r - (long) PPM_GETR(*pP) ) * (FS_SCALE/16);
		    thisrerr[col + 2] += err * 7;
		    nextrerr[col    ] += err * 3;
		    nextrerr[col + 1] += err * 5;
		    nextrerr[col + 2] += err;
		    err = ( g - (long) PPM_GETG(*pP) ) * (FS_SCALE/16);
		    thisgerr[col + 2] += err * 7;
		    nextgerr[col    ] += err * 3;
		    nextgerr[col + 1] += err * 5;
		    nextgerr[col + 2] += err;
		    err = ( b - (long) PPM_GETB(*pP) ) * (FS_SCALE/16);
		    thisberr[col + 2] += err * 7;
		    nextberr[col    ] += err * 3;
		    nextberr[col + 1] += err * 5;
		    nextberr[col + 2] += err;
		    }
		else
		    {
		    err = ( r - (long) PPM_GETR(*pP) ) * (FS_SCALE/16);
		    thisrerr[col    ] += err * 7;
		    nextrerr[col + 2] += err * 3;
		    nextrerr[col + 1] += err * 5;
		    nextrerr[col    ] += err;
		    err = ( g - (long) PPM_GETG(*pP) ) * (FS_SCALE/16);
		    thisgerr[col    ] += err * 7;
		    nextgerr[col + 2] += err * 3;
		    nextgerr[col + 1] += err * 5;
		    nextgerr[col    ] += err;
		    err = ( b - (long) PPM_GETB(*pP) ) * (FS_SCALE/16);
		    thisberr[col    ] += err * 7;
		    nextberr[col + 2] += err * 3;
		    nextberr[col + 1] += err * 5;
		    nextberr[col    ] += err;
		    }
		}

	    if ( ( ! floyd ) || fs_direction )
		{
		++col;
		++pP;
		}
	    else
		{
		--col;
		--pP;
		}
	    }
	while ( col != limitcol );

	if ( floyd )
	    {
	    register long *temperr;

	    temperr = thisrerr;
	    thisrerr = nextrerr;
	    nextrerr = temperr;

	    temperr = thisgerr;
	    thisgerr = nextgerr;
	    nextgerr = temperr;

	    temperr = thisberr;
	    thisberr = nextberr;
	    nextberr = temperr;

	    fs_direction = ! fs_direction;
	    }

	ppm_writeppmrow( stdout, pixelrow, cols, maxval, 0 );
	}
    }


/* Nice color-searching routine by James H. Bruner. */
static int
find_best_color( pixel* pP, pixval maxval, colorhist_vector chv, int newcolors )
    {
    register int i;
    register int gdist;
    register int r, g, b;
    register int gidx, low, mid, high;
    long dist, newdist;
    int idx;
    static long* squares = (long*) 0;
    static long* squares_o;
    static int* rs;
    static int* gs;
    static int* bs;

    if ( squares == (long*) 0 )
	{
	/* Initialize.  Cache squares of [-maxval..maxval]. */
	squares = (long*) pm_allocrow( 2 * (int) maxval + 1, sizeof(long) );
	squares_o = squares + maxval;
	for ( i = 0; i <= (int) maxval; ++i )
	    squares_o[i] = squares_o[-i] = (long) i * (long) i;

	/* Sort by the green component.  This is just a guess but since
	** green contributes most to luminosity, I think that there might be
	** better distribution across green.  It really doesn't matter which
	** color though.
	*/
	qsort(
	    (char*) chv, newcolors, sizeof(struct colorhist_item),
	    green_compare );

	/* Cache colormap values. */
	rs = (int*) pm_allocrow( newcolors, sizeof(int) );
	gs = (int*) pm_allocrow( newcolors, sizeof(int) );
	bs = (int*) pm_allocrow( newcolors, sizeof(int) );
	for ( i = 0; i < newcolors; ++i )
	    {
	    rs[i] = PPM_GETR( chv[i].color );
	    gs[i] = PPM_GETG( chv[i].color );
	    bs[i] = PPM_GETB( chv[i].color );
	    }
	}

    r = PPM_GETR( *pP );
    g = PPM_GETG( *pP );
    b = PPM_GETB( *pP );

    /* The colormap is sorted by the green component.  Do binary search. */
    low = 0;
    high = newcolors - 1;
    while ( low < high )
	{
	mid = ( high + low ) / 2;
	gdist = gs[mid];
	if ( gdist == g )
	    high = low = mid;
	else
	    if ( gdist < g )
		low = mid + 1;
	    else
		high = mid - 1;
	}
    gidx = high;
    if ( gidx < 0 )
	gidx = 0;
	/* Could be <0 if value less than 1st.  Can't be
	** greater than newcolors because division for mid
	** rounds down.
	*/

    /* The color map is sorted by the green component.  We've searched and now
    ** have the index at the closest green match.  This means, that we can now
    ** search up and down from this location and only have to search until
    ** the green component alone is worse than the best total color match.
    ** Because it is sorted by green, we know it will only get worse!
    */

    /* First search down. */
    dist = 2000000000;
    for ( i = gidx; i >= 0; --i )
	{
	gdist = squares_o[g - gs[i]];
	if ( gdist >= dist )
	    break;
	newdist = gdist + squares_o[r - rs[i]] + squares_o[b - bs[i]];
	if ( newdist < dist )
	    {
	    dist = newdist;
	    idx = i;
	    }
	}

    /* Now search up. */
    for ( i = gidx + 1; i < newcolors; ++i )
	{
	gdist = squares_o[g - gs[i]];
	if ( gdist >= dist )
	    break;
	newdist = gdist + squares_o[r - rs[i]] + squares_o[b - bs[i]];
	if ( newdist < dist )
	    {
	    dist = newdist;
	    idx = i;
	    }
	}

    return idx;
    }