*Anton works for Truda Software, a Fortran and C consulting firm. He can be reached through the DDJ offices.*

Note that the RGB components of a color can be viewed as lying along the axes of a *color space*, or RGB cube; see Figure 1. For a 24-bit image, the color space is (practically speaking) continuous, since the smallest difference between colors is imperceptible. This continuous color space must be mapped to 256 discrete colors, which are then used to represent the color space. Mapping a continuous variable to a discrete set of values is called *quantization*.

A simple approach to color quantization is to take the 256 most-common colors in the true-color image. For each color in the input image, we then search for the closest color in the set of 256 colors, and display that color instead. Actually, the 256 colors are loaded in the display system's color map or video look-up table (LUT), and the pixels' LUT indexes are written out to the display hardware. The algorithm that picks the 256 most-common colors to load into the LUT is the *popularity algorithm*. Its simplicity comes at a price, because it gives poor results for many images. Other color-quantization algorithms may give excellent results, but may be very time consuming. For example, as Wan et al. report, it "may take more than 20 hours on a VAX 780 computer to produce 256 clusters for a full-color image" using an algorithm known as the *K-means algorithm*. Color quantization is a post-processing procedure on an image, and should normally take only a fraction of the time to generate the image. Thus, such run times are in most cases unacceptable.

Some algorithms quantize an image without regard to its content; their dominant feature is speed. They use a fixed number of standard colors. For example, 256 colors require eight bits of color per pixel, and the bits are divided among red, green, and blue. Since the human eye is less sensitive to blue, red and green each get three bits, and blue two bits. This gives eight shades of red, eight of green, and four of blue. Other colors are mixtures of these primaries. The problem is that for a particular image, most standard colors may never be used. For example, if an image contains a lot of red and very little blue, all the entries assigned to blue are effectively wasted. This is a fundamental flaw of image-independent quantizers--if the image and the quantizer don't match, the results are poor.

For example, Table 1 is a histogram for 14 pixels that have six unique colors; for clarity, the colors have no blue component. The task is to apply the median-cut algorithm to this histogram and find four output colors.

The initial rectangle is in Figure 2(a); crosses are used to indicate the pixels' colors. The distance from maximum to minimum along the *R*-axis is 80--5=75; along the *G*-axis, 80--20=60. Thus, the longest distance from maximum to minimum is along the *R*-axis. The median along this axis is *R*=30, so this is where the rectangle is split, and two tight-fitting rectangles are placed around the two new regions; see Figure 2(b).

When applying the procedure to the left rectangle, the axis with the longest distance is the *G*-axis. The median along this axis is *G*=50--this is where the axis is split. For the right rectangle, the axis with the longest distance from maximum to minimum is also the *G*-axis. The median along this axis is *G*=40, and this is where the axis is split. The final rectangles are shown in Figure 2(c), where each rectangle has about the same number of pixels. Following this, the centroids of the colors in each rectangle are computed. The resulting values are the output colors; see Figure 2(d).

There are two approaches to remapping an input image: *fast remapping* and *best remapping*. With fast remaps, the centroid of a cube represents all the colors enclosed by the cube. This is often good enough, but does not give the best results because the centroid of the cube in which a color falls may not be the closest centroid. Best remap searches the list of output colors for the closest color, but because searching is involved, it's slower.

The median-cut algorithm is normally described as recursive. However, it's easier and more practical to implement it as a nonrecursive routine. (Any recursive algorithm can be converted to an iterative algorithm.) To see why, apply the recursive algorithm in Figure 3(a), assuming it's used to generate four cubes.

The first division of the RGB cube results in two smaller cubes called *CubeA _{1}* and

Another problem an implementation should address is the special case in which a cube encloses only one color, typically repeated many times. Obviously, such a cube cannot be split. One method of dealing with this is to split one of the other cubes; otherwise, you end up with one less cube than you want. To do this, you need a list of all the current cubes. However, this isn't available with the algorithms in Figure 3(a) and Figure 3(c), since you implicitly let the algorithm save the cubes on the program stack during the recursive calls. You can deal with this by increasing the maximum allowable level by 1 when encountering a cube with a single color. Unfortunately, this again opens the door for the problem of dividing smaller cubes while a larger cube is still waiting on the stack. Furthermore, the situation where the desired number of final cubes is not a power of 2, further complicates the process.

To address these problems, I used the following nonrecursive method. A list of the cubes computed so far is maintained in an array. During each iteration, the list of cubes is scanned for the cube with the smallest level, but cubes with single colors are ignored. Several cubes are normally candidates for splitting, and you could devise a rule for picking the best one to split, but I simply used the first one found. This cube is then split at the median, perpendicular to the longest axis. This increases the number of cubes by one. One of the new cubes is saved in the slot of the cube just split, and the other is added to the bottom of the list of cubes. When the desired number of cubes is generated, or all of the cubes enclose a single color, the splitting phase terminates. This algorithm is summarized in Figure 3(d).

At this point you have a list of cubes, and the next step is to compute their centroids. These colors are the output colors, or output color map. The histogram is no longer required, and you can use its space for a look-up table that holds the output colors--an *inverse color map*. The inverse color map works just like the histogram, where 24-bit colors are mapped to 15-bit colors that serve as indexes into the histogram. However, the inverse color map contains the input colors' indexes in the color map, instead of their count.

I've seen implementations of the median-cut algorithm that compute the inverse color map for the whole RGB cube. However, this is unnecessary because you know which colors each cube encloses, and they're the only ones required. Initialization of the inverse color map depends on the kind of remapping desired. With fast remap, initialization is accomplished by replacing the counts for all the colors in a cube with the centroid of the cube. With best remap, initialization is accomplished by replacing the counts for all the colors in a cube by the *closest* centroid, found by searching the list of centroids.

I've used the C standard-library function *qsort* (based on the Quicksort algorithm with its average run-time complexity of *O(N log _{2}N)*) to do the sort. However, the algorithm has a potentially serious pitfall. When it's called to sort data that is already in order, the run-time complexity is O(N

Another problem with Quicksort is that it is often implemented as a recursive routine, and can rapidly deplete the program stack when a large number of data points must be sorted. A nonrecursive implementation needs much less auxiliary storage. Thus, a good implementation of Quicksort is normally nonrecursive, and it randomizes the input data somewhat to provide for the case were the data is already in sorted order. Many implementations, however, don't do this. For example, I examined the Microsoft C 5.1 run-time library source, and its *qsort* is a fairly simple recursive implementation of Quicksort. For a production version of the median-cut algorithm, you might want to replace the *qsort* routine with a sort routine that's less efficient on the average, but has a better worst-case performance. Heapsort has both an average and worst-case run-time complexity of *O(N log _{2}N)*, and it needs no (or very little) auxiliary storage.

The histogram, which I call *Hist*, is accessed indirectly via the array *HistPtr* that functions as an index into the histogram. *HistPtr* contains the position of the colors in the actual histogram. When the histogram (or parts of it) must be sorted, *HistPtr* (or parts of it) is sorted. A structure is used to represent a cube. In Listing One, this structure (*cube_t*) has several members, but the essential members are *lower *and *upper, *and they point to the opposite corners of a cube in *HistPtr*.

To see how this works, return to the data in Table 1, and rework the example in terms of the *Hist* and *HistPtr*. Also, keep an eye on Figure 2 to see the correspondence between the data structures and the geometric interpretation of the median-cut algorithm. Assume that after the histogram is constructed, it looks like *Hist*; see Figure 4(a). Note the "holes" in the histogram--this is typical. Also, this example assumes that the colors are initially sorted on green. C_{1} has the smallest green component, C_{4} has the second smallest, and so on. Just after *Hist* is constructed, *HistPtr* is filled as in Figure 4(a). The initial cube is the whole RGB cube, so that *Cube.lower* is set to 0, and *Cube.upper* to 5. This completes the initialization.

Now the cube must be split. The red axis is the longest, so we sort the colors along this axis. To do this, sort the array *HistPtr* instead of *Hist*. To find the median along the axis, start at *Cube.lower*; its value is 0. We look in *HistPtr[0]*; its value is 5, and this corresponds to C_{2} in *Hist*. C_{2}'s count is 4, and this is the running sum. Now we look at *Cube.lower*+1's corresponding color, C_{0}, and the running sum becomes 3+4=7. This is half the total for the cube, so it is split to form two cubes, *CubeA _{1}*, and

Next we split *CubeB _{1}*. The colors in

To compute the centroid for, say, *CubeB _{2}*, start at

MedianCut, the program in Listing One was developed on and for DOS PCs. To ease porting, I've tried to adhere to ANSI C. I also used *typedef*s for some variables when the number of bits was important. When moving to a different platform, you might want to change these so that they reflect the sizes of the target system. The code compiles and executes cleanly with Microsoft's 5.1 compiler (large memory model). If the /W4 compiler switch is used with version 6.0 of this compiler, some warning messages are issued, but they can safely be ignored.

MedianCut takes three arguments: the image histogram *Hist*, which contains the 15-bit colors; *maxcubes*, the desired number of output colors (this can be any number between 1 and 256, but the upper limit can be changed by altering the *#define MAXCOLORS* preprocessor directive); and a *maxcubes *x 3 array *ColMap* that MedianCut fills with the output colors, where the red component of color *i* is *r=ColMap[i][0]*; that of green, *g=ColMap[i][1]; *and that of *blue*, *b=ColMap[i][2]*. MedianCut returns the number of actual colors, which may be less than *maxcubes *if the input image already contains less colors than the requested number.

The *#define FAST_REMAP* preprocessor directive controls whether the inverse color map is initialized with the fast-remap or the best-remap method. If this directive is deleted or commented out, initialization is done according to the best-remap method.

Additionally, I've written a typical driver routine for MedianCut that, in this case, converts a Targa type-2 image file to a Targa type 1. Truevision's Targa type-2 image format is a popular format for true-color images, and Targa type-1 images are color-mapped, so to convert from type 1 to type 2, we must quantize the true-color image. This program is available electronically; see page 3.

Instead of dynamically allocating the space for *HistPtr* and the array for the list of cubes, I've chosen to have these static variables. These arrays occupy about 70 Kbytes. The histogram occupies another 64 Kbytes, and depending on your compiler's *qsort* routine, several Kbytes of stack may be required. Add to this about another 20 Kbytes for the file buffers, as well as the space required for the program code, and a total of about 200 Kbytes of memory are required.

Figures 5 and 6 each show two images before and after application of the median-cut algorithm (fast remap). Both were generated with the POV-Ray (Persistence of Vision) ray-tracing program (see "Ray Tracing and the POV-Ray Toolkit," by Craig A. Lindley, *DDJ*, July 1994). The 256-color version of "grid" is indistinguishable from the true-color original, which shows how effective the median-cut algorithm can be. However, with the "lamp" image, there are some false contours in areas of low contrast, but this may not show up in the photographs, since it is quite difficult to reproduce a video display faithfully in print.

This implementation of the median-cut algorithm is quite fast. With the fast-remap method, it takes about eight seconds to quantize 640x480 true-color versions of the images in Figures 5 and 6 on a 20-MHz AT clone. Much of this time is I/O, and the actual color quantization takes less than three seconds.

Heckbert, P. "Color Image Quantization for Frame Buffer Display." *Computer Graphics* (July 1982).

Figure 1 The RGB cube.

Figure 2 The median-cut algorithm: (a) original rectangle; (b) after splitting once; (c) after splitting three times; (d) output colors.

(a) Split(Cube){ if (ncubes == 4) return; find longest axis of Cube; cut Cube at median to form CubeA, CubeB; Split(CubeA); Split(CubeB); } (b) Figure (c) maxlevel = 2; Split(Cube,level){ if (ncubes == 4) return; if (Cube's level == maxlevel) return; find longest axis of Cube; cut Cube at median to form CubeA, CubeB; Split(CubeA, level+1); Split(CubeB, level+1); } (d) build initial cube from histogram; set initial cube's level to 0; insert initial cube in list of cubes; ncubes = 1; while (ncubes < maxcubes){ search for Cube with smallest level; find the longest axis of Cube; find the median along this axis; cut Cube at median to form CubeA, CubeB; set CubeA's level = Cube's level + 1; set CubeB's level = Cube's level + 1; insert CubeA in Cube's slot; add CubeB to end of list of cubes; ncubes = ncubes + 1; }

Figure 4 (a) Hist and HistPtr after initialization; (b) sorted on red and split at the median; (c) CubeB_{1}, sorted on green; (d) CubeB_{1}, split at the median.

Figure 5 Grid: (a) original, true-color image; (b) 256-color version (image-file courtesy of Dan Farmer).

Figure 6 Lamp: (a) original, true-color image; (b) 256-color version (anonymous image-file description).

Color (r,g)-coordinates Count C_{0}(20,40) 3 C_{1}(40,20) 2 C_{2}(5,60) 4 C_{3}(50,80) 2 C_{4}(60,30) 1 C_{5}(80,50) 2

Cube HistPtr.lower HistPtr.upper Colors Enclosed Cube Centroid A_{3}0 0 C_{0}(20,40) B_{3}1 1 C_{2}(5,60) A_{2}2 3 C_{1},C_{4}(46.7,23.3) B_{2}4 5 C_{5},C_{3}(65,65)

Example 1 Computing cube centroids.

/* median.c -- Anton Kruger, Copyright (c) Truda Software, 215 Marengo Rd, ** #2, Oxford, IA 52322-9383 ** Description: Contains an implementation of Heckbert's median-cut color ** quantization algorithm. ** Compilers: MSC 5.1, 6.0. ** Note: 1) Compile in large memory model. 2) Delete "#define FAST_REMAP" ** statement below in order to deactivate fast remapping. */ #define FAST_REMAP #include <stdio.h> #include <stddef.h> /* for NULL */ #include <stdlib.h> /* for "qsort" */ #include <float.h> /* for FLT_MAX, FLT_MIN */ #define MAXCOLORS 256 /* maximum # of output colors */ #define HSIZE 32768 /* size of image histogram */ typedef unsigned char byte; /* range 0-255 */ typedef unsigned short word; /* range 0-65,535 */ typedef unsigned long dword; /* range 0-4,294,967,295 */ /* Macros for converting between (r,g,b)-colors and 15-bit */ /* colors follow. */ #define RGB(r,g,b) (word)(((b)&~7)<<7)|(((g)&~7)<<2)|((r)>>3) #define RED(x) (byte)(((x)&31)<<3) #define GREEN(x) (byte)((((x)>>5)&255)<< 3) #define BLUE(x) (byte)((((x)>>10)&255)<< 3) typedef struct { /* structure for a cube in color space */ word lower; /* one corner's index in histogram */ word upper; /* another corner's index in histogram */ dword count; /* cube's histogram count */ int level; /* cube's level */ byte rmin,rmax; byte gmin,gmax; byte bmin,bmax; } cube_t; static cube_t list[MAXCOLORS]; /* list of cubes */ static int longdim; /* longest dimension of cube */ static word HistPtr[HSIZE]; /* points to colors in "Hist" */ void Shrink(cube_t * Cube); void InvMap(word * Hist, byte ColMap[][3],word ncubes); int compare(const void * a1, const void * a2); word MedianCut(word Hist[],byte ColMap[][3], int maxcubes) { /* Accepts "Hist", a 32,768-element array that contains 15-bit color counts ** of input image. Uses Heckbert's median-cut algorithm to divide color ** space into "maxcubes" cubes, and returns centroid (average value) of each ** cube in ColMap. Hist is also updated so that it functions as an inverse ** color map. MedianCut returns the actual number of cubes, which may be ** less than "maxcubes". */ byte lr,lg,lb; word i,median,color; dword count; int k,level,ncubes,splitpos; void *base; size_t num,width; cube_t Cube,CubeA,CubeB; /* Create the initial cube, which is the whole RGB-cube. */ ncubes = 0; Cube.count = 0; for (i=0,color=0;i<=HSIZE-1;i++){ if (Hist[i] != 0){ HistPtr[color++] = i; Cube.count = Cube.count + Hist[i]; } } Cube.lower = 0; Cube.upper = color-1; Cube.level = 0; Shrink(&Cube); list[ncubes++] = Cube; /* Main loop follows. Search the list of cubes for next cube to split, which ** is the lowest level cube. A special case is when a cube has only one ** color, so that it cannot be split. */ while (ncubes < maxcubes){ level = 255; splitpos = -1; for (k=0;k<=ncubes-1;k++){ if (list[k].lower == list[k].upper) ; /* single color */ else if (list[k].level < level){ level = list[k].level; splitpos = k; } } if (splitpos == -1) /* no more cubes to split */ break; /* Must split the cube "splitpos" in list of cubes. Next, find longest ** dimension of cube, and update external variable "longdim" which is ** used by sort routine so that it knows along which axis to sort. */ Cube = list[splitpos]; lr = Cube.rmax - Cube.rmin; lg = Cube.gmax - Cube.gmin; lb = Cube.bmax - Cube.bmin; if (lr >= lg && lr >= lb) longdim = 0; if (lg >= lr && lg >= lb) longdim = 1; if (lb >= lr && lb >= lg) longdim = 2; /* Sort along "longdim". This prepares for the next step, namely finding ** median. Use standard lib's "qsort". */ base = (void *)&HistPtr[Cube.lower]; num = (size_t)(Cube.upper - Cube.lower + 1); width = (size_t)sizeof(HistPtr[0]); qsort(base,num,width,compare); /* Find median by scanning through cube, computing a running sum. When ** running sum equals half the total for cube, median has been found. */ count = 0; for (i=Cube.lower;i<=Cube.upper-1;i++){ if (count >= Cube.count/2) break; color = HistPtr[i]; count = count + Hist[color]; } median = i; /* Now split "Cube" at median. Then add two new cubes to list of cubes.*/ CubeA = Cube; CubeA.upper = median-1; CubeA.count = count; CubeA.level = Cube.level + 1; Shrink(&CubeA); list[splitpos] = CubeA; /* add in old slot */ CubeB = Cube; CubeB.lower = median; CubeB.count = Cube.count - count; CubeB.level = Cube.level + 1; Shrink(&CubeB); list[ncubes++] = CubeB; /* add in new slot */ if ((ncubes % 10) == 0) fprintf(stderr,"."); /* pacifier */ } /* We have enough cubes, or we have split all we can. Now compute the color ** map, inverse color map, and return number of colors in color map. */ InvMap(Hist, ColMap,ncubes); return((word)ncubes); } void Shrink(cube_t * Cube) { /* Encloses "Cube" with a tight-fitting cube by updating (rmin,gmin,bmin) ** and (rmax,gmax,bmax) members of "Cube". */ byte r,g,b; word i,color; Cube->rmin = 255; Cube->rmax = 0; Cube->gmin = 255; Cube->gmax = 0; Cube->bmin = 255; Cube->bmax = 0; for (i=Cube->lower;i<=Cube->upper;i++){ color = HistPtr[i]; r = RED(color); if (r > Cube->rmax) Cube->rmax = r; if (r < Cube->rmin) Cube->rmin = r; g = GREEN(color); if (g > Cube->gmax) Cube->gmax = g; if (g < Cube->gmin) Cube->gmin = g; b = BLUE(color); if (b > Cube->bmax) Cube->bmax = b; if (b < Cube->bmin) Cube->bmin = b; } } void InvMap(word * Hist, byte ColMap[][3],word ncubes) { /* For each cube in list of cubes, computes centroid (average value) of ** colors enclosed by that cube, and loads centroids in the color map. Next ** loads histogram with indices into the color map. A preprocessor directive ** #define FAST_REMAP controls whether cube centroids become output color ** for all the colors in a cube, or whether a "best remap" is followed. */ byte r,g,b; word i,j,k,index,color; float rsum,gsum,bsum; float dr,dg,db,d,dmin; cube_t Cube; for (k=0;k<=ncubes-1;k++){ Cube = list[k]; rsum = gsum = bsum = (float)0.0; for (i=Cube.lower;i<=Cube.upper;i++){ color = HistPtr[i]; r = RED(color); rsum += (float)r*(float)Hist[color]; g = GREEN(color); gsum += (float)g*(float)Hist[color]; b = BLUE(color); bsum += (float)b*(float)Hist[color]; } /* Update the color map */ ColMap[k][0] = (byte)(rsum/(float)Cube.count); ColMap[k][1] = (byte)(gsum/(float)Cube.count); ColMap[k][2] = (byte)(bsum/(float)Cube.count); } #ifdef FAST_REMAP /* Fast remap: for each color in each cube, load the corresponding slot ** in "Hist" with the centroid of the cube. */ for (k=0;k<=ncubes-1;k++){ Cube = list[k]; for (i=Cube.lower;i<=Cube.upper;i++){ color = HistPtr[i]; Hist[color] = k; } if ((k % 10) == 0) fprintf(stderr,"."); /* pacifier */ } #else /* Best remap: for each color in each cube, find entry in ColMap that has ** smallest Euclidian distance from color. Record this in "Hist". */ for (k=0;k<=ncubes-1;k++){ Cube = list[k]; for (i=Cube.lower;i<=Cube.upper;i++){ color = HistPtr[i]; r = RED(color); g = GREEN(color); b = BLUE(color); /* Search for closest entry in "ColMap" */ dmin = (float)FLT_MAX; for (j=0;j<=ncubes-1;j++){ dr = (float)ColMap[j][0] - (float)r; dg = (float)ColMap[j][1] - (float)g; db = (float)ColMap[j][2] - (float)b; d = dr*dr + dg*dg + db*db; if (d == (float)0.0){ index = j; break; } else if (d < dmin){ dmin = d; index = j; } } Hist[color] = index; } if ((k % 10) == 0) fprintf(stderr,"."); /* pacifier */ } #endif return; } int compare(const void * a1, const void * a2) { /* Called by the sort routine in "MedianCut". Compares two ** colors based on the external variable "longdim". */ word color1,color2; byte C_{1},C_{2}; color1 = (word)*(word *)a1; color2 = (word)*(word *)a2; switch (longdim){ case 0: C_{1}= RED(color1), C_{2}= RED(color2); break; case 1: C_{1}= GREEN(color1), C_{2}= GREEN(color2); break; case 2: C_{1}= BLUE(color2), C_{2}= BLUE(color2); break; } return ((int)(C_{1}-C_{2})); }