/* 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" #define pm_closer(x) pm_close(x) #define pm_closew(x) pm_close(x) #define pm_srandom(x) srandom(x) #define pm_random(x) random(x) /* This define controls how big the histogram can get before the program ** rescales and tries again. */ /* #define MAXCOLORS 32767 */ #define MAXCOLORS 50000 /* 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 compact 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 compact 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 compact 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. */ 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 ARGS(( colorhist_vector chv, int colors, long sum, pixval maxval, int newcolors )); static void choose_reps ARGS(( box_vector bv, int newcolors, colorhist_vector chv, pixval maxval )); static void bound_box ARGS(( struct box* boxP, colorhist_vector chv )); static void add_box_to_array ARGS(( pixel**** color_array, struct box* box, pixval maxval )); static void add_pixel_to_array ARGS(( pixel**** color_array, int neighbor_dist, int r, int g, int b, pixel* newpP, pixval maxval )); static void map_colors ARGS(( 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 )); static int find_best_color ARGS(( pixel* pP, pixval maxval, colorhist_vector chv, int newcolors )); static int red_compare ARGS(( colorhist_vector ch1, colorhist_vector ch2 )); static int green_compare ARGS(( colorhist_vector ch1, colorhist_vector ch2 )); static int blue_compare ARGS(( colorhist_vector ch1, colorhist_vector ch2 )); void main( argc, argv ) int argc; char* argv[]; { FILE* ifp; pixel** pixels; pixel** mappixels; pixel* pixelrow; int readlines; int argn, rows, cols, format, maprows, mapcols, i; pixval mapmaxval, maxval; int newcolors, colors; colorhist_vector chv; int floyd, writemap; pixel**** color_array; char* usage = "[-floyd|-fs] [ppmfile]\n [-floyd|-fs] -map mapfile [ppmfile]"; ppm_init( &argc, argv ); argn = 1; floyd = 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", 2 ) || pm_keymatch( argv[argn], "-nofloyd", 2 ) ) floyd = 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 histogram, making the new colormap. */ pm_message( "choosing %d color%s...", newcolors, (newcolors == 1) ? "" : "s" ); bv = median_cut( chv, colors, rows * cols, maxval, newcolors ); ppm_freecolorhist( chv ); /* Choose representative colors for the boxes. */ choose_reps( bv, newcolors, chv, maxval ); /* 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 compact 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 ** compact 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 compact 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 colormap into a 3-D compact 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 ); } } /* 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 ); /* All done. */ pm_closew( stdout ); exit( 0 ); } /* 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. */ #if __STDC__ static box_vector median_cut( colorhist_vector chv, int colors, long sum, pixval maxval, int newcolors ) #else /*__STDC__*/ static box_vector median_cut( chv, colors, sum, maxval, newcolors ) colorhist_vector chv; int colors, newcolors; long sum; pixval maxval; #endif /*__STDC__*/ { 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. */ #if __STDC__ static void choose_reps( box_vector bv, int newcolors, colorhist_vector chv, pixval maxval ) #else /*__STDC__*/ static void choose_reps( bv, newcolors, chv, maxval ) box_vector bv; int newcolors; colorhist_vector chv; pixval maxval; #endif /*__STDC__*/ { 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 compact array of pixel*, indexed ** by [r][g][b]. ** ** This is used to convert a set of boxes produced by median-cut into ** the initial compact array. */ #if __STDC__ static void add_box_to_array( pixel**** color_array, struct box* box, pixval maxval ) #else /*__STDC__*/ static void add_box_to_array( color_array, box, maxval ) pixel**** color_array; struct box* box; pixval maxval; #endif /*__STDC__*/ { 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 compact array of pixel*, indexed ** by [r][g][b]. ** ** This is used both to convert a user-supplied colormap into the initial ** compact array, and during Floyd-Steinberg error propagation to add ** slightly modified entries to the array. */ #if __STDC__ static void add_pixel_to_array( pixel**** color_array, int neighbor_dist, int r, int g, int b, pixel* newpP, pixval maxval ) #else /*__STDC__*/ static void add_pixel_to_array( color_array, neighbor_dist, r, g, b, newpP, maxval ) pixel**** color_array; int neighbor_dist; int r, g, b; pixel* newpP; pixval maxval; #endif /*__STDC__*/ { 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 #if __STDC__ 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 ) #else /*__STDC__*/ static void map_colors( pixels, pixelrow, readlines, ifp, rows, cols, maxval, format, color_array, chv, newcolors, floyd ) 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; #endif /*__STDC__*/ { long *thisrerr; long *nextrerr; long *thisgerr; long *nextgerr; long *thisberr; long *nextberr; int row, col; int fs_direction; int limitcol; 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; } 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 += thisrerr[col + 1] / FS_SCALE; g += thisgerr[col + 1] / FS_SCALE; b += thisberr[col + 1] / FS_SCALE; if ( r < 0 ) r = 0; else if ( r > maxval ) r = maxval; if ( g < 0 ) g = 0; else if ( g > maxval ) g = maxval; if ( b < 0 ) b = 0; else if ( b > maxval ) b = maxval; PPM_ASSIGN( *pP, r, g, b ); } /* Check closest-color compact 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. */ 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. */ #if __STDC__ static int find_best_color( pixel* pP, pixval maxval, colorhist_vector chv, int newcolors ) #else /*__STDC__*/ static int find_best_color( pP, maxval, chv, newcolors ) pixel* pP; pixval maxval; colorhist_vector chv; int newcolors; #endif /*__STDC__*/ { 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; } static int red_compare( ch1, ch2 ) colorhist_vector ch1, ch2; { return (int) PPM_GETR( ch1->color ) - (int) PPM_GETR( ch2->color ); } static int green_compare( ch1, ch2 ) colorhist_vector ch1, ch2; { return (int) PPM_GETG( ch1->color ) - (int) PPM_GETG( ch2->color ); } static int blue_compare( ch1, ch2 ) colorhist_vector ch1, ch2; { return (int) PPM_GETB( ch1->color ) - (int) PPM_GETB( ch2->color ); }