watson@ames.arpa (John S. Watson) (12/19/87)
Hi Folks,
I'm looking for either references to, or code for a really
nice color dithering algorithm. The algorithm produces 8-bit dithered
images that are almost indistinguishable from the original 24-bit images.
The algorithm colors a pixels by somehow taking into account the
quantization errors in coloring the neighboring pixels.
I think a similar algorithm is used in Sun's NeWS.
I think the algorithm was first described by Floyd and Steinberg
in 1975 but I don't have the exact reference(s). If I had some harder
references I could probably code up the algorithm myself, but I also
don't want to re-invent the wheel.
--
John S. Watson ARPA: watson@ames.arpa
NASA Ames Research Center UUCP: ...!ames!watson
Any opinions expressed herein are solely the responsibility of the
author and do not represent the opinions of NASA or the U.S. Governmentlou@bearcat.rutgers.edu (Lou Steinberg) (12/22/87)
In article <3703@ames.arpa> watson@ames.arpa (John S. Watson) writes: > I'm looking for either references to, or code for a really > nice color dithering algorithm. The algorithm produces 8-bit dithered > images that are almost indistinguishable from the original 24-bit images. > > The algorithm colors a pixels by somehow taking into account the > quantization errors in coloring the neighboring pixels. > I think a similar algorithm is used in Sun's NeWS. > > I think the algorithm was first described by Floyd and Steinberg > in 1975 but I don't have the exact reference(s). If I had some harder > references I could probably code up the algorithm myself, but I also > don't want to re-invent the wheel. I think you are referring to the article Paul Heckert, "Color Image Quantization for Frame Buffer Display", Computer Graphics, V. 16 #3, July 1982, p. 297. The dithering method he uses is an adaptation of an algorithm developed by Bob Floyd and myself for B/W halftoning. Our algorithm is described in Robert Floyd and Louis Steinberg, "An Adaptive Algorithm for Spatial Greyscale", Proceedings of the Society for Information Display, V.17 #2, 2nd Quarter 1976, P. 75 A good recent source on the B/W problem is Digital Halftoning, a book by Robert Ulichney, MIT Press, 1987 (based on his PhD thesis). Among other things, it discusses some improvements to our algorithm. The basic idea is that reducing a 24 bit color to an 8 bit color introduces some error. The closest available 8 bit color might have a certain amount too much red, a certain amount too little green, and a certain amount too much blue (or whatever other dimensions you are using for your color space). The algorithm scans the image. At each pixel it chooses the closest 8 bit color and calculates the error it has thus introduced. Some fraction of this error is used to adjust the color of each of those neighboring pixels that have not been processed yet. E.g. if the 8 bit color is too red, reduce the red value of the neighboring pixels. Heckert's article also discusses how to choose the specific 256 colors to use for the 8 bit picture (he assumes you have an 8 bit -> 24 bit color map in your display), and how to efficiently find the closest color in this set to any given 24 bit color. I believe Heckert is at Pixar now. Perhaps he can provide the code, although it should not be hard to rewrite it yourself. Also, Paul Roetling at Xerox's labs in Webster (Mass??) has been working on similar things. I can provide the code for the B/W version, with the standard disclaimer that "the secretary will disavow any knowledge ...". It is in C, runs on Suns and should run on anything. I do not know what is used in NeWS, but the B/W images I've seen do look as if they were produced by our algorithm. I'd be curious to know what method they use, since it seems to run very fast. -- Lou Steinberg uucp: {pretty much any major site}!rutgers!aramis.rutgers.edu!lou arpa: lou@aramis.rutgers.edu
lou@bearcat.rutgers.edu (Lou Steinberg) (12/22/87)
Oops - in my previous message I misspelled Paul Heckbert's last name.
My appologies, Paul.
--
Lou Steinberg
uucp: {pretty much any major site}!rutgers!aramis.rutgers.edu!lou
arpa: lou@aramis.rutgers.educfchiesa@bsu-cs.UUCP (Christopher F. Chiesa) (12/30/87)
In article <545@bearcat.rutgers.edu>, lou@bearcat.rutgers.edu (Lou Steinberg) writes: [... mucho del stuffo deleted...] > > I believe Heckert is at Pixar now. Perhaps he can provide the code, > although it should not be hard to rewrite it yourself. Also, Paul > Roetling at Xerox's labs in Webster (Mass??) has been working on ^^^^^^^^^^^^^^^^ Xerox Webster is in Webster, NY - a suburb of Rochester. Chris Chiesa Senior, CS Dept., Ball State University, Muncie IN (Native of Rochester, NY) !rutgers!iuvax!bsu-cs!cfchiesa
lee@rocksanne.UUCP (Lee Moore) (01/05/88)
> > I believe Heckert is at Pixar now. Perhaps he can provide the code, > > although it should not be hard to rewrite it yourself. Also, Paul > > Roetling at Xerox's labs in Webster (Mass??) has been working on > ^^^^^^^^^^^^^^^^ > Xerox Webster is in Webster, NY - a suburb of Rochester. > In case anybody is interested, Paul is reachable electronically as roetling.wbst128@xerox.com [] [] [] [] [] -- Lee Moore -- Xerox Webster Research Center, birthplace of the XGP UUCP: {seismo, allegra, decvax, cmcl2, topaz}!rochester!rocksanne!lee Arpa Internet: Moore.wbst@xerox.arpa DDN: +1 (716) 422-2496
garyo@masscomp.UUCP (Gary Oberbrunner) (01/07/88)
In article <3703@ames.arpa> watson@ames.UUCP (John S. Watson) writes: > I'm looking for either references to, or code for a really >nice color dithering algorithm. The algorithm produces 8-bit dithered >images that are almost indistinguishable from the original 24-bit images. Floyd & Steinberg's dithering algorithm is described in the following paper: Floyd, R.W. and Steinberg, L. "An Adaptive Algorithm for Spatial Gray Scale." SID 75m Int. Symp. Dig. Tech. Papers (1975), p. 36. I don't actually have this paper, just a reference to it in Heckbert's "Color Image Quantization for Frame Buffer Display (Siggraph proceedings, 1982, p. 297). You probably want to check this paper out too. But the way F/S dither works is like this (it's so simple...) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - First you come up with a color map - use uniform segmentation for speed or Heckbert's method for quality. Then you work from the top left to the bottom right of the image (pixel order), as follows: Find the closest match in the color map to the true pixel color. Use that color map entry to represent that pixel. Compute the r,g,b error resulting from the above approximation. Add 3/8 of that error to the pixel to the right, and 3/8 to the pixel below the current pixel. The remaining 1/4 goes to the pixel diagonally below. Don't forget to clip the rgb, as well as check your boundary conditions. And that's it - compute the next pixel! No dither maps, no muss no fuss no bother. And it's real fast since you can do the 3/8 and 1/4 by shifting. On some images it can help to add some random noise; for instance if you have a smooth (flat-shaded) surface that's just below one of the colors in the map, the error will build up slowly until you get a line of brighter color somewhere in the surface. A bit o' noise gets rid of that line nicely. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Of course choosing the color map is the fun part. :-) :-) Gary Oberbrunner Here is a hack C program which implements both Heckbert's algorithm and F/S dither (no claims of any kind are made for this program - in fact this article ended above with my name, and this is all line nois)@(#*&... No seriously, this code needs real help - it was written for a one-shot deal, and it is very inefficient in many places. But it does work ok. The worst part is where every possible color mapping is computed - SLOW! And Remap() ain't a speed-demon either, but that's where Floyd et. al. come in... Enjoy (watch out for possible signature at end...) ---------------------------------------------------------------------------- Remember, -Truth is not beauty; Information is not knowledge; / Beauty is not love; Gary Oberbrunner Knowledge is not wisdom; / Love is not music; ...!masscomp!garyo Wisdom is not truth; ----/ Music is the best. - FZ ---------------------------------------------------------------------------- ---------------------- c u t h e r e ---------------------------------- #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create: # quantize.c (most everything in here) # split.c (the color-space-splitting algorithm) # quantize.h (structs & externs) # std.h (my "gotta have it" stuff) # This archive created: Wed Jan 6 21:58:18 1988 export PATH; PATH=/bin:/usr/bin:$PATH if test -f 'quantize.c' then echo shar: "will not over-write existing file 'quantize.c'" else cat << \SHAR_EOF > 'quantize.c' /*---------------------------------------------------------------------- * Color image quantizer, from Paul Heckbert's paper in * Computer Graphics, vol.16 #3, July 1982 (Siggraph proceedings), * pp. 297-304. * By Gary Oberbrunner, copyright c. 1988. *---------------------------------------------------------------------- */ #include <stdio.h> #include <fcntl.h> #include <std.h> #include "quantize.h" static char SCCSId[] = "SCCS ID: %W% %G% %U"; int NCOLS = 0; char *progname; int debug = 1; int dither = 0; /* This turns on dithering */ int fs_dither = 1; /* This turns on Floyd-Steinberg dithering */ #define DITH 8 int Dither[DITH][DITH] = { { 0, 32, 8, 40, 2, 34, 10, 42}, {48, 16, 56, 24, 50, 18, 58, 26}, {12, 44, 4, 36, 14, 46, 6, 38}, {60, 28, 52, 20, 62, 30, 54, 22}, { 3, 35, 11, 43, 1, 33, 9, 41}, {51, 19, 59, 27, 49, 17, 57, 25}, {15, 47, 7, 39, 13, 45, 5, 37}, {63, 31, 55, 23, 61, 29, 53, 21} }; /* This is only temporarily global: */ out_pix_t *quant_image; /* the quantized image - cmap indices */ main (argc, argv) int argc; char **argv; { char *filename, *outfilename; int rows = 0, NCOLS = 0, cmaplen = 0; rgb_pix_t *image; /* rgb image storage */ cmap_t *cmap; /* color map; rgb byte triples */ unsigned int *hist; /* histogram of color values: indices are formed from 5 bits each of r,g,b */ unsigned int *maptable; /* quantized color map table: indices are colors, values are color map indices */ int i; /* temp! */ progname = argv[0]; if (argc != 6) { dbug1("Usage: %s infile rows cols colormaplen outfile\n", argv[0]); dbug0("Converts 24-bit rgb image to n-bit colormapped image.\n"); exit(1); } filename = argv[1]; rows = atoi(argv[2]); NCOLS = atoi(argv[3]); cmaplen = atoi(argv[4]); outfilename = argv[5]; /* dynamically allocate all the various structures */ NEW(image, rgb_pix_t, rows * NCOLS * RGB_BYTES); NEW(quant_image, out_pix_t, rows * NCOLS); NEW(cmap, cmap_t, 3 * cmaplen); NEW(hist, unsigned int, HIST_SIZE); maptable = hist; /* reuse the same storage for the map table */ /* Perform the algorithm */ ReadImage(filename, image, rows, NCOLS); MakeHist(image, rows * NCOLS, hist); /* ShowHist(hist);*/ MakeCmap(hist, cmap, cmaplen); ComputeMapping(cmap, cmaplen, maptable); RemapImage(image, rows, NCOLS, cmap, cmaplen, maptable, quant_image); WriteImage(outfilename, quant_image, rows, NCOLS, cmap, cmaplen); /* Show the result */ #ifdef GRAPHICS mgiasngp(0,0); for (i = 0; i < cmaplen*3; i+=3) { printf("Before LoadMap: cmap[%4d] = %3d, %3d, %3d\n", i/3, cmap[i], cmap[i+1], cmap[i+2]); } LoadMap(cmap, cmaplen); Draw(quant_image, rows, NCOLS); mgideagp(); #endif } /* This assumes that the image is stored in packed r,g,b byte triples in */ /* row-major order. */ ReadImage(name, array, rows, cols) char *name; rgb_pix_t *array; int rows, cols; { int file; int xfercount; file = open(name, O_RDONLY); if (file == -1) { dbug2("%s: Can't open file %s for reading.\n", progname, name); perror(name); exit(1); } dbug0("Reading file..."); xfercount = read(file, array, RGB_BYTES * rows * cols); if (xfercount == -1) { dbug2("%s: Can't read file %s.\n", progname, name); perror(name); exit(1); } if (xfercount != RGB_BYTES * rows * cols) { dbug2("File %s ended prematurely after %d bytes: continuing.\n", name, xfercount); } dbug0("Done.\n"); } /* Write out the quantized image */ WriteImage(name, array, rows, cols, cmap, cmaplen) char *name; out_pix_t *array; int rows, cols; cmap_t *cmap; int cmaplen; { int file; int xfercount; struct FileData { int rows, cols, cmaplen, dummy; } fdata; int i; fdata.rows = rows; fdata.cols = cols; fdata.cmaplen = cmaplen; fdata.dummy = 0; file = open(name, O_WRONLY | O_TRUNC | O_CREAT, 0666); if (file == -1) { dbug2("%s: Can't open file %s for writing.\n", progname, name); perror(name); exit(1); } dbug0("Writing file..."); dbug1("Header is %d bytes.\n", sizeof(fdata)); xfercount = write(file, &fdata, sizeof(fdata)); if (xfercount != sizeof(fdata)) { dbug3("%s: Can't write file %s's header. Wrote %d bytes.\n", progname, name, xfercount); if (xfercount == -1) perror(name); exit(1); } dbug1("Color map is %d bytes.\n", 3 * cmaplen); xfercount = write(file, cmap, 3 * cmaplen); if (xfercount != 3*cmaplen) { dbug3("%s: Can't write file %s's color map. Wrote %d bytes.\n", progname, name, xfercount); if (xfercount == -1) perror(name); exit(1); } xfercount = write(file, array, rows * cols); if (xfercount != rows * cols) { dbug3("%s: Can't write image to file %s. Wrote %d bytes.\n", progname, name, xfercount); if (xfercount == -1) perror(name); exit(1); } dbug0("Done.\n"); close(file); } MakeHist(image, npix, hist) rgb_pix_t *image; int npix; unsigned int *hist; { register int pix; /* pixel index */ register rgb_pix_t *ppix; /* pointer to the actual triple */ /* Clear out the histogram array */ for (pix = HIST_SIZE - 1; pix >= 0; pix--) hist[pix] = 0; dbug0("Making histogram...\n"); for (pix = npix-1; pix >= 0; pix--) { if (pix % 25000 == 0) dbug1("%d \r", pix); ppix = PIXPIX(image, pix); /* Increment the histogram bucket corresponding to the value at pix: Only five bits of the color value are used to increase clumping in the histogram. The histogram is indexed by the five bits each of red, green and blue all strung together. */ hist[(*ppix >> HIST_SHIFT << (HIST_QUANT_BITS * 2)) | (*(ppix+1) >> HIST_SHIFT << HIST_QUANT_BITS) | (*(ppix+2) >> HIST_SHIFT)]++; } dbug0("Done.\n"); } ShowHist(hist) unsigned int *hist; { register int i, nmaps = 0; printf("RED\tGREEN\tBLUE\n"); for (i = 0; i < HIST_SIZE; i++) { if (hist[i] != 0) { nmaps++; printf("%-d\t%-d\t%-d\t= %d.\n", i >> 10 << 3, (i >> 5 & 0x1f) << 3, (i & 0x1f) << 3, hist[i]); } } dbug1("%d different pixel values (5-bit quantized).\n", nmaps); } /* Compute the map color for each original color */ /* This should be done dynamically while remapping... */ ComputeMapping(cmap, cmaplen, maptable) cmap_t *cmap; int cmaplen; unsigned int *maptable; { register int mapcolor, d, mindist, minmap; int orig_color; for (orig_color = 0; orig_color < HIST_SIZE; orig_color++) { mindist = BIG_DISTANCE; for (mapcolor = 3*(cmaplen-1); mapcolor >= 0; mapcolor -= 3) { if ((d = DIST(((orig_color >> 10) & 0x1f) << 3, ((orig_color >> 5) & 0x1f) << 3, (orig_color & 0x1f) << 3, cmap[mapcolor], cmap[mapcolor+1], cmap[mapcolor+2])) < mindist) { mindist = d; minmap = mapcolor; } } maptable[orig_color] = minmap/3; if ((orig_color & 0x3ff) == 0) dbug1("%d/31 \r", orig_color >> 10); } } /* Pick color map values that are close to the real ones, using the map * computed above in ComputeMapping(). */ /* Now with optional dither! */ RemapImage(image, rows, cols, cmap, cmaplen, maptable, quant_image) rgb_pix_t *image; /* original image (rgb) */ int rows, cols; /* size of the image */ cmap_t *cmap; /* color map (rgb) */ int cmaplen; /* number of entries in cmap */ unsigned int *maptable; /* mapping table from rgb5 to cmap indices */ out_pix_t *quant_image; /* quantized image */ { int npix = rows * cols; register unsigned int pix; register int mapentry; int nearestcm; dbug0("Remapping image..."); for (pix = 0; pix < npix; pix++) { mapentry = ((PIXRED(image, pix) >> 3 & 0x1f) << 10) | ((PIXGREEN(image, pix) >> 3 & 0x1f) << 5) | (PIXBLUE(image, pix) >> 3 & 0x1f); ASSERT (mapentry >= 0 && mapentry < HIST_SIZE) nearestcm = maptable[mapentry]; if (fs_dither) { register int r, g, b; register int errr, errg, errb; /* Compute error in current pixel */ errr = PIXRED(image, pix) - cmap[nearestcm * 3]; errg = PIXGREEN(image, pix) - cmap[nearestcm * 3 + 1]; errb = PIXBLUE(image, pix) - cmap[nearestcm * 3 + 2]; /* Propagate errors right and down */ if (pix % cols != cols - 1) { /* right gets 3/8 of the error */ r = (int) PIXRED(image,pix+1) + errr * 3 / 8; g = (int) PIXGREEN(image,pix+1) + errg * 3 / 8; b = (int) PIXBLUE(image,pix+1) + errb * 3 / 8; PIXRED(image,pix+1) = clamp(r, 0, 255); PIXGREEN(image,pix+1) = clamp(g, 0, 255); PIXBLUE(image,pix+1) = clamp(b, 0, 255); } if (pix < npix - cols) { /* down gets another 3/8 */ r = (int) PIXRED(image,pix+cols) + errr * 3 / 8; g = (int) PIXGREEN(image,pix+cols) + errg * 3 / 8; b = (int) PIXBLUE(image,pix+cols) + errb * 3 / 8; PIXRED(image,pix+cols) = clamp(r, 0, 255); PIXGREEN(image,pix+cols) = clamp(g, 0, 255); PIXBLUE(image,pix+cols) = clamp(b, 0, 255); } if (pix < npix - cols - 1) { /* and diagonally down gets the last 1/4. */ r = (int) PIXRED(image,pix+cols+1) + errr / 4; g = (int) PIXGREEN(image,pix+cols+1) + errg / 4; b = (int) PIXBLUE(image,pix+cols+1) + errb / 4; PIXRED(image,pix+cols+1) = clamp(r, 0, 255); PIXGREEN(image,pix+cols+1) = clamp(g, 0, 255); PIXBLUE(image,pix+cols+1) = clamp(b, 0, 255); } /* Store the current pixel */ quant_image[pix] = nearestcm; } else if (dither) { register unsigned int r, g, b; register unsigned int nr, ng, nb; register int oppr, oppg, oppb; int oppnrcm; /* true rgb value for pixel */ r = PIXRED(image, pix); g = PIXGREEN(image, pix); b = PIXBLUE(image, pix); /* rgb values for nearest color */ nr = cmap[nearestcm * 3]; ng = cmap[nearestcm * 3 + 1]; nb = cmap[nearestcm * 3 + 2]; /* Color as far from rgb as nrngnb but in the other direction */ oppr = min(255, max(0, 2 * (int) r - (int) nr)); oppg = min(255, max(0, 2 * (int) g - (int) ng)); oppb = min(255, max(0, 2 * (int) b - (int) nb)); /* Nearest match for opposite color: */ oppnrcm = maptable[((oppr >> 3 & 0x1f) << 10) | ((oppg >> 3 & 0x1f) << 5) | (oppb >> 3 & 0x1f)]; /* If they're not the same, dither between them. */ /* Dither constant is measured by where the true color lies between the two nearest approximations. Since the most nearly opposite color is not necessarily on the line from the nearest through the true color, some triangulation error can be introduced. In the worst case the r-nr distance can actually be less than the nr-oppr distance. */ if (oppnrcm != nearestcm) { register int dither_const; int oppcmr = cmap[oppnrcm * 3]; int oppcmg = cmap[oppnrcm * 3 + 1]; int oppcmb = cmap[oppnrcm * 3 + 2]; dither_const = min(63, (64 * DIST(r, g, b, nr, ng, nb)) / DIST(nr, ng, nb, oppcmr, oppcmg, oppcmb)); if (!(dither_const >= 0 && dither_const < 64)) { dbug2("Pix: %d, Dither constant: %d\n", pix, dither_const); fprintf(stderr,"Near: %d %d %d\tTrue: %d %d %d\tOpp: %d %d %d / %d %d %d\n", nr, ng, nb, r, g, b, oppr, oppg, oppb, oppcmr, oppcmg, oppcmb); } ASSERT (dither_const >= 0 && dither_const < 64); if (Dither[(pix/cols)%DITH][(pix%cols)%DITH] < dither_const) quant_image[pix] = oppnrcm; else quant_image[pix] = nearestcm; } else quant_image[pix] = nearestcm; } else { quant_image[pix] = nearestcm; } } dbug0("Done.\n"); /* Cmap is still ok here */ } /* Load up the frame buffer color map - device specific */ LoadMap(cmap, cmaplen) cmap_t *cmap; int cmaplen; { register int i; dbug1("Loading %d color map entries...", cmaplen); #ifdef GRAPHICS for (i = 0; i < cmaplen*3; i+=3) { register unsigned int r = cmap[i], g = cmap[i+1], b = cmap[i+2]; mgicm(i/3, (r << 16) | (g << 8) | b); dbug2("Color %d is 0x%x.\n", i/3, (r << 16) | (g << 8) | b); } #endif dbug0("done.\n"); } /* Draw the remapped image - device specific */ Draw(image, rows, cols) out_pix_t *image; int rows, cols; { #ifdef GRAPHICS mgiimage(image, 0, 0, cols-1, rows-1); #endif } SHAR_EOF fi if test -f 'split.c' then echo shar: "will not over-write existing file 'split.c'" else cat << \SHAR_EOF > 'split.c' /*---------------------------------------------------------------------- * split up color space for quantize program *---------------------------------------------------------------------- */ #include <stdio.h> #include <std.h> #include "quantize.h" /* Pick your choice of algorithm: */ /* MIDPOINT_CUT or POPULARITY */ /* (Choose in makefile) */ #ifdef MIDPOINT_CUT typedef unsigned int color_index_t; /* My method here is to allocate a clist array big enough for all the * possible colors, and then recursively split it. Since the recursion * is depth-first and we're building a LIST, not a TREE of boxes, we can * scramble the colors within a parent box because that parent won't be used * any more. We end up with a list of boxes into which all the existing * colors are partitioned, as described in Heckbert's paper. */ struct Box { struct Box *next; int rmin, rmax, gmin, gmax, bmin, bmax; color_index_t *clist; /* an array of color indices */ int clist_len; /* how many cindexes in the list */ }; static int nboxes, maxboxes; /* number of boxes so far and max needed */ MakeCmap(hist, cmap, cmaplen) unsigned int *hist; cmap_t *cmap; int cmaplen; { register int i; struct Box *box; nboxes = 1; maxboxes = cmaplen; NEW(box, struct Box, 1); dbug0("Making color boxes.\n"); MakeTopBox(box, hist); /* Make the first (top-level) box */ while (nboxes < maxboxes) SplitBoxes(box); FindColors(box, cmap, cmaplen); } MakeTopBox(box, hist) /* Make the top box from the histogram */ struct Box *box; color_index_t *hist; { register int i; ASSERT(box != NULL); ASSERT(hist != NULL); /* Allocate enough color cells to hold all the colors */ NEW(box->clist, color_index_t, HIST_SIZE); box->clist_len = 0; box->next = NULL; /* For each possible rgb color */ for (i = 0; i < HIST_SIZE; i++) { if (hist[i] > 0) /* if any pixels are that color, */ box->clist[box->clist_len++] = i; } Compact(box); fprintf(stderr,"Top box at 0x%x: %d colors, R (%d -> %d) G (%d -> %d) B (%d -> %d)\n", box, box->clist_len, box->rmin, box->rmax, box->gmin, box->gmax, box->bmin, box->bmax); } /* Shrink the box limits until it just encloses the points within it. */ Compact(box) struct Box *box; { register int i; ASSERT(box != NULL); box->rmin = 256; box->gmin = 256; box->bmin = 256; /* Greater than any value */ box->rmax = -1; box->gmax = -1; box->bmax = -1; /* Less than any rgb value */ for (i = box->clist_len-1; i >= 0; i--) { int r = (box->clist[i] >> 10) & 0x1f; /* colors associated with 'i' */ int g = (box->clist[i] >> 5) & 0x1f; int b = (box->clist[i]) & 0x1f; /* Expand the box until it just fits the colors tightly */ if (r < box->rmin) box->rmin = r; if (g < box->gmin) box->gmin = g; if (b < box->bmin) box->bmin = b; if (r > box->rmax) box->rmax = r; if (g > box->gmax) box->gmax = g; if (b > box->bmax) box->bmax = b; } } /* Split all the boxes once each, stop if we've got enough */ SplitBoxes(box) struct Box *box; { struct Box *nextbox; dbug1("Splitting all the boxes with %d done so far.\n", nboxes); while (box != NULL && nboxes < maxboxes) { nextbox = box->next; SplitBox(box); /* SplitBox changes box->next! */ box = nextbox; } } static char longdim; /* char representing the longest dimension */ int CompareColors(c1, c2) color_index_t *c1, *c2; { switch (longdim) { case 'r': return ((*c1 >> 10 & 0x1f) - (*c2 >> 10 & 0x1f)); case 'g': return ((*c1 >> 5 & 0x1f) - (*c2 >> 5 & 0x1f)); case 'b': return ((*c1 & 0x1f) - (*c2 & 0x1f)); default: fprintf(stderr,"Illegal longdim 0x%x in CompareColors.\n", longdim); exit(1); } } /* Divide the box across its longest dimension at its median point * This is NOT recursive - we split all the boxes once, then if there's * more colors to be assigned we split all of them again, and so on. */ SplitBox(box) struct Box *box; { struct Box *newbox; /* the 'other half' of the split box */ ASSERT(box != NULL); if (nboxes > maxboxes) return; if (box->clist_len <= 1) return; nboxes++; if (box->rmax - box->rmin >= box->gmax - box->gmin && box->rmax - box->rmin >= box->bmax - box->bmin) longdim = 'r'; else if (box->gmax - box->gmin >= box->rmax - box->rmin && box->gmax - box->gmin >= box->bmax - box->bmin) longdim = 'g'; else if (box->bmax - box->bmin >= box->gmax - box->gmin && box->bmax - box->bmin >= box->rmax - box->rmin) longdim = 'b'; else { fprintf(stderr,"No longest dimension in SplitBox??!\n"); fprintf(stderr,"This box: R (%d -> %d) G (%d -> %d) B (%d -> %d)\n", box->rmin, box->rmax, box->gmin, box->gmax, box->bmin, box->bmax); exit(1); } qsort(box->clist, box->clist_len, sizeof(color_index_t), CompareColors); /* Set up the new box and put it after the old box */ NEW(newbox, struct Box, 1); newbox->next = box->next; box->next = newbox; /* If len is odd, put the extra one in the old box */ newbox->clist_len = box->clist_len / 2; box->clist_len -= newbox->clist_len; newbox->clist = box->clist + box->clist_len; Compact(box); Compact(newbox); } /* Find the representative color in each box and put it into cmap. */ FindColors(box, cmap, cmaplen) /* here box is a list of boxes */ struct Box *box; cmap_t *cmap; int cmaplen; { register int r, g, b; register int i; int cml = 0; while (box != NULL) { r = 0, g = 0, b = 0; for (i=0; i < box->clist_len; i++) { r += (box->clist[i] >> 10 & 0x1f) << 3; g += (box->clist[i] >> 5 & 0x1f) << 3; b += (box->clist[i] & 0x1f) << 3; } ASSERT(r/box->clist_len < 256); ASSERT(g/box->clist_len < 256); ASSERT(b/box->clist_len < 256); cmap[cml*3] = r / box->clist_len; cmap[cml*3+1] = g / box->clist_len; cmap[cml*3+2] = b / box->clist_len; box = box->next; cml++; } ASSERT(cmaplen == cml); /* since there should be 'cmaplen' boxes */ for (i = 0; i < cmaplen*3; i+=3) { printf("Cmap[%4d] = %3d, %3d, %3d\n", i/3, cmap[i], cmap[i+1], cmap[i+2]); } } #endif #ifdef POPULARITY static unsigned int *HiSt; CompareHist(v1, v2) int *v1, *v2; { return (HiSt[*v2] - HiSt[*v1]); } MakeCmap(hist, cmap, cmaplen) unsigned int *hist; cmap_t *cmap; int cmaplen; { register int i; unsigned int *hptrs; /* This is the popularity algorithm (Boyle & Lippman, Archmac, 1978) */ /* First we sort the histogram: set up a bunch of pointers and sort them */ dbug1("Getting %d ints for hptrs...", HIST_SIZE); /* NEW(hptrs, unsigned int, HIST_SIZE);*/ hptrs = (unsigned int *)quant_image; /* use existing storage temporarily */ dbug0("initting it..."); for (i = 0; i < HIST_SIZE; i++) hptrs[i] = i; HiSt = hist; dbug0("Sorting..."); qsort (hptrs, HIST_SIZE, sizeof(unsigned int), CompareHist); dbug0("making color map..."); for (i = 0; i < cmaplen; i++) { cmap[3*i] = hptrs[i] >> 10 << 3; cmap[3*i+1] = (hptrs[i] >> 5 & 0x1f) << 3; cmap[3*i+2] = (hptrs[i] & 0x1f) << 3; } dbug0("done.\n"); dbug1("Max pixels in one bucket: %d.\n", hist[hptrs[0]]); } #endif POPULARITY SHAR_EOF fi if test -f 'quantize.h' then echo shar: "will not over-write existing file 'quantize.h'" else cat << \SHAR_EOF > 'quantize.h' /*---------------------------------------------------------------------- * header file for color-quantizer * by Gary Oberbrunner 1/5/88 *---------------------------------------------------------------------- */ #define HIST_QUANT_BITS 5 /* Histogram clump quantization */ #define HIST_SHIFT (8 - HIST_QUANT_BITS) #define HIST_SIZE (1 << (HIST_QUANT_BITS * 3)) #define RGB_BYTES 4 /* Bytes per rgb pixel (4 for alpha channel) */ /* We need the next definition to reference a variable-size array */ /* NCOLS is the number of columns in the array */ /* Each pixel is stored as an rgb byte triplet */ /* This is the address of the pixel triplet: */ #define RCPIX(array, row, col) (array + RGB_BYTES * (row * NCOLS + col)) #define PIXPIX(array, pixel) (array + RGB_BYTES * (pixel)) /* These are the pixel components: */ #define RED(array, row, col) (*(RCPIX(array, row, col))) #define GREEN(array, row, col) (*(RCPIX(array, row, col) + 1)) #define BLUE(array, row, col) (*(RCPIX(array, row, col) + 2)) #define PIXRED(array, pix) (*(PIXPIX(array, pix))) #define PIXGREEN(array, pix) (*(PIXPIX(array, pix) + 1)) #define PIXBLUE(array, pix) (*(PIXPIX(array, pix) + 2)) /* Here is the rgb-space distance metric: */ #define DIST(r1,g1,b1,r2,g2,b2) \ (int) \ (3 * ((r1)-(r2)) * ((r1)-(r2)) + \ 4 * ((g1)-(g2)) * ((g1)-(g2)) + \ 2 * ((b1)-(b2)) * ((b1)-(b2))) #define BIG_DISTANCE 1000000 /* bigger than any real distance */ typedef unsigned char cmap_t; typedef unsigned char rgb_pix_t; typedef unsigned char out_pix_t; extern out_pix_t *quant_image; /* the quantized image - cmap indices */ extern int debug; SHAR_EOF fi if test -f 'std.h' then echo shar: "will not over-write existing file 'std.h'" else cat << \SHAR_EOF > 'std.h' /* Standard include file with lots of generally useful stuff */ /* by Gary Oberbrunner */ #define min(x,y) ((x) < (y) ? (x) : (y)) #define max(x,y) ((x) > (y) ? (x) : (y)) #define clamp(x, mn, mx) (((x) <= (mn)) ? (mn) : (((x) >= (mx)) ? (mx) : (x))) #define SECURITY #ifdef SECURITY #define ASSERT(exp) { if (!(exp)) { fprintf(stderr,\ "Assertion error: %s, line %d. Assertion exp failed.\n",\ __FILE__, __LINE__);\ exit(69); } } #else #define ASSERT(exp) #endif #define dbug0(str) \ if (debug) fprintf(stderr, str); #define dbug1(str, arg1) \ if (debug) fprintf(stderr, str, arg1); #define dbug2(str, arg1, arg2) \ if (debug) fprintf(stderr, str, arg1, arg2); #define dbug3(str, arg1, arg2, arg3) \ if (debug) fprintf(stderr, str, arg1, arg2, arg3); #define dbug4(str, arg1, arg2, arg3, arg4) \ if (debug) fprintf(stderr, str, arg1, arg2, arg3, arg4); extern char *malloc(); /* * NEW is a macro to malloc 'n' new variables of type 'type'. */ #define NEW(var, type, num) \ if ((var = (type *) malloc((num) * sizeof(type)))==NULL) \ { fprintf(stderr,\ "ERROR: Out of memory.\nLast request:\n");\ fprintf(stderr,\ "\tNumber: num\n\tType : type\n\tName : var.\n");\ fprintf(stderr,"File %s, line %d\n", __FILE__, __LINE__);\ fprintf(stderr,\ "\tTotal requested: num=%d x sizeof(type)=%d\n\t =%d bytes.\n",\ num, sizeof(type), (num) * sizeof(type));\ exit(99); } SHAR_EOF fi exit 0 # End of shell archive -- Remember, -Truth is not beauty; Information is not knowledge; / Beauty is not love; Gary Oberbrunner Knowledge is not wisdom; / Love is not music; ...!masscomp!garyo Wisdom is not truth; ----/ Music is the best. - FZ
awpaeth@watcgl.waterloo.edu (Alan W. Paeth) (01/09/88)
A comment on dither - the "quick and dirty" approach to uniform division of
the color space (so that R, G and B can be treated separably) very often
slices up 8 bit pixels as 3 bits for red and green, 2 bits for blue (where
the eye is least sensitive).
This is an unnecessary oversimplification that leaves blue with only two
mid-range intensities. It suggests itself because we tend to think in terms
of "bits" per color at the hardware level, not "discrete number of intensity
steps". This approach is further suggested because the pixel channels for
R G and B can be written separably through the proper setting of the write
mask, but in practice, the pixel is normally written as a single byte anyway,
not as three consecutive "inking" passes. One useful property of this approach
is that there are a number of "greys" in the color space, as the channel
precisions are all multiples of each other.
A better approach is to form a color table containing a Cartesian product of
intensites on a range that is not necessarily a power of two. For example, if
we take red = [0..5], green = [0..6], blue = [0..4] (suitably normalized to
represent intensities on the range [0.0..1.0] then there are 6*7*5 = 210 color
table entries, and noticibly better representation in blue.
Our Graphics Lab takes this approach for writing on our VAX GPX's under X.
We must use a uniform space, because many simultaneous image (windows) might
be on the display, and each would have a different histogram (if Heckbert's
approach were used), potentially overflowing the color table. In our case, the
X-write IM tool takes color images of arbitrary precision and Floyd-Steinberg
converts them into this space of 6,7,5 steps, and thus all color images share
identical entries within the color table (which is a consequence of how X
manages requests for color table entries). B/W images typically allocate a 32
level grey ramp which can co-exist within the remaining space, leaving
14 (256=210+32+14) entries for use by other applications, should they require
colors other than saturated yellow, blue, white, etc.
On GPX's with smaller color tables, the same approach is used. Once again, ~3/4
of the hardware color table is allocated, giving 12 free color entries, which
the program divides as 2*3*2 in RGB. Other useful sets have been 6*6*6, which
gives the largest uniform allocation in 256 entries, thus allowing for exact
greys, and 6*7*6, which uses 252 out of 256 colors and is much better overall
than the 8*8*4 set which is used implicitly when one allocates using the old
"three bits red and green, two for blue" scheme.
/Alan Paeth
Computer Graphics Lab
University of Waterloolou@aramis.rutgers.edu (Lou Steinberg) (01/11/88)
In article <2741@masscomp.UUCP> garyo@masscomp.UUCP (Gary Oberbrunner) writes: > [for] Floyd & Steinberg's dithering algorithm [...] > On some images it can help to add some random noise; for instance if you > have a smooth (flat-shaded) surface that's just below one of the colors in > the map, the error will build up slowly until you get a line of brighter > color somewhere in the surface. A bit o' noise gets rid of that line nicely. The other way to get rid of it is to alternate scanning left to right and right to left, i.e. even scan lines one way, odd scan lines the other. The line of brighter color essentially comes because pixels on that line do not communicate - there is no path for the error value at one to affect the value chosen for the other. By scanning in a zig zag fashion, you ensure that each pixel communicates with every other pixel scanned after it. -- Lou Steinberg uucp: {pretty much any major site}!rutgers!aramis.rutgers.edu!lou arpa: lou@aramis.rutgers.edu
lou@aramis.rutgers.edu (Lou Steinberg) (01/11/88)
In article <2741@masscomp.UUCP> garyo@masscomp.UUCP (Gary Oberbrunner) writes: > But the wayd [Floyd/Steinberg] dither works is like this (it's so simple...) [...] > Add 3/8 of that error to the pixel to the right, and 3/8 to the pixel > below the current pixel. The remaining 1/4 goes to the pixel diagonally > below. I should point out that the actual fractions we used were, assuming you are at X, moving left to right: X 7/16 3/16 5/16 1/16 Note that the error goes to four neighbors, not three. I think this will probably do better (at least for black and white) than the 3/8-3/8-1/4 distribution, at the cost of greater processing. I have seen the 3/8-3/8-1/4 distribution described as "our" algorithm before, but I have no idea who the credit really belongs to. Also, I should add that if you do zig-zag scanning (see my immediately previous message), it is sufficient (but not quite as good) to send half the error one pixel ahead (e.g. to the right on lines you scan left to right), and half one pixel straight down. Again, this is for black and white; I've not tried it with color. -- Lou Steinberg uucp: {pretty much any major site}!rutgers!aramis.rutgers.edu!lou arpa: lou@aramis.rutgers.edu
mmt@dciem.UUCP (Martin Taylor) (01/22/88)
>A comment on dither - the "quick and dirty" approach to uniform division of >the color space (so that R, G and B can be treated separably) very often >slices up 8 bit pixels as 3 bits for red and green, 2 bits for blue (where >the eye is least sensitive). > >This is an unnecessary oversimplification that leaves blue with only two >mid-range intensities. >[There follows a discussion of coding colour tables] It is not strictly true that the eye is least sensitive to colour variation in the neighbourhood of blue. In fact, the MacAdam discrimninatio ellipses are smallest near blue and largest near green. What IS true is that the eye does not discriminate blue-yellow contrasts in small areas, nor does it have any blue-yellow edge detectors. Hence, blue is a good candidate for colour dither. I have never tried this, having not been associated with colour display since about 12 years ago, but the following seems like a good idea, based on visual function: Give as much colour resolution per pixel as you can to red-green contrast, and minimize dither in that part of the colour space. For blue, which dominates the blue-yellow contrast of most regions, use either 1 or 2 bits, but average over substantial regions by dithering. If one dithers over a 4X4 region with 1 bit, one has 16 blue levels over regions that probably are not far from the resolution limit of the eye for blue-yellow contrast. Unfortunately, this will not be the whole story, because the blue phosphor contributes more to brightness than does a good central blue, but it should be a reasonable approximation. As a practical matter, I suspect that the traditional 3+3+2 red-green-blue distribution of bits is not bad if one dithers red and green on a 3-pixel region, and blue on a 9-pixel region. The 3-pixel regions would be L-shaped or inverted-L-shaped, interlocking for either red or green, with the interlock column offset 1-pixel between red and green. At discrete edges in the image, don't dither, but use the best per-pixel value. The existence of the edge will mask any subtle colour error. These are untested ideas, but off-hand I can't see why they would not give good subtle shading in a colour space of 9x9x18 (red-blue-green) over gradually varying regions of a picture, plus good edge resolution in other regions. You actually want better blue colour resolution and worse blue spatial resolution than for red and green, which is what this scheme would give you. But I have no idea how computationally expensive it would be to have these non-rectangular dither areas which had different boundaries for all three primary colours. If anyone tries this, I would like to know how it works (and to get credit if it is magnificent!). -- Martin Taylor {allegra,linus,ihnp4,floyd,ubc-vision}!utzoo!dciem!mmt {uw-beaver,qucis,watmath}!utcsri!dciem!mmt mmt@zorac.arpa Magic is just advanced technology ... so is intelligence. Before computers, the ability to do arithmetic was proof of intelligence. What proves intelligence now?
awpaeth@watcgl.waterloo.edu (Alan W. Paeth) (01/26/88)
The discussion was a follow-on to my previous posting, in which I suggest splitting 256 color map entries as a Cartesian product of RGB in which the Green had the greatest precision, because of the eye's sensitivity to green. >It is not strictly true that the eye is least sensitive to colour >variation in the neighbourhood of blue. In fact, the MacAdam discrimninatio >ellipses are smallest near blue and largest near green. What IS true >is that the eye does not discriminate blue-yellow contrasts in small >areas, nor does it have any blue-yellow edge detectors. Hence, blue >is a good candidate for colour dither... (discussion of some possible dithering techniques) >...You actually want better blue colour resolution >and worse blue spatial resolution than for red and green, which is >what this scheme would give you. But I have no idea how computationally >expensive it would be to have these non-rectangular dither areas which >had different boundaries for all three primary colours. In reverse order: [2] This really hits the nail on the head. The preliminary design of the "Dandelion" by Butler Lampson at Xerox PARC in the late 70's went for a video chain which had r3g3b2 precision, but the b2 channel was clocked in such a fashion so as to produce pixels of four bits resolution (in pixel quantization), but of double the width, thus losing some spatial resolution. This makes sense, as the eye not only has better contrast discernment on the red-green channel (as was pointed out), but cannot even focus well on blue, being from -1D to -2D out as the wavelength decreases, owing to chromatic aberration. Such a viedo design wouldn't be that hard a retrofit for many systems, as the total number of blue bits clocked per scanline (and thus, per screen) remains unchanged -- just a two bit buffer is needed to build up the full blue pixel from the memory "subpixels", plus some logic gates. [1] Yes, the MacAdam ellipses (and those by Stiles) demonstrate empirically that hue changes are most discernable for BRG, in that order, even when the perceptual hue non-linearities of the CIE x'y' chart (on which the ellipses are often drawn) are taken into account. When dithering a subject of largely varying blue, we choose this order of bit allocation -- as for instance when doing a sunset against a blue sky. In general, though, the high blue precision is most relevant only when a strictly blue color is to be rendered; the extra bits might not allow for a large perceptual color shift when the other two primaries are active. The real problem here stems from treating the creation of a suitable color space. The previous posting suggests an opponent space. I've tried Floyd- Steinberg techniques in this space, as well as LUV, YIQ and HSV with varying (and noticeably distinct) results. These spaces are largely a change in basis, or straight-forward non-linear transformation of RGB. As such, they are highly symmetric, so forcing a degree of high precision to some corner of the space by bit allocation almost always gives rise to a "conjugate" space of increased precision in an area which is uninteresting (the precision increase in this conjugate area is not very perceptual and is thus wasted). What is really needed is a Floyd-Steinburg algorithm which dithers against a candidate list of target output colors. These colors would be defined to be either perceptually "equidistant", or alternately they would be derived from the input image by using a more general histogramming technique (cluster analysis), which would not treat the input primaries as separate channels. Although this potentially means histogramming against 2^24 bins, note that most images seldom exceed 2^20 pixels in size, a reduction of 16:1. The real difficulty with such a program is doing the "nearest colored pixel" function, given, say, an input of r8g8b8 data and an output of 256 distinct pixels. For small output (eg, 16 colors for dithering onto a PC), one can search the output pixel space in linear fashion. For output pixels of large precision, say, 16 bits, a r5g6b5 separable Floyd scheme with equalization (Heckbert, SIGGRAPH '80) is more than adequate. The interesting case is when the output device is of intermediate size -- in practice this is also the most common case, as when the output device has an upper limit of 256 output pixel colors. This size is too large for linear brute force, but small enough that sublinear algorithms with asymptotically better bounds (probably involving 3D Voroni diagrams) are large and not necessarily faster given the expense in building and probing the colorspace datastructure. They are also a royal pain to code. /Alan Paeth University of Waterloo Computer Graphics Laboratory