Color quantization

Color quantization is the process of reducing number of colors used in an image while trying to maintain the visual appearance of the original image. In general, it is a form of cluster analysis, if each RGB color value is considered as a coordinate triple in the 3D colorspace. There are some well know algorithms [1], each with its own advantages and drawbacks.

Task: Take an RGB color image and reduce its colors to some smaller number (< 256). For this task, use the frog as input and reduce colors to 16, and output the resulting colors. The chosen colors should be adaptive to the input image, meaning you should not use a fixed palette such as Web colors or Windows system palette. Dithering is not required.

Note: the funny color bar on top of the frog image is intentional.

C

Using an octree to store colors. Here are only the relevant parts. For full C code see Color_quantization/C. It's different from the standard octree method in that:

Each node can both be leaf node and have child nodes;
Leaf nodes are not folded until all pixels are in. This removes the possibility of early pixels completely bias the tree. And child nodes are reduced one at a time instead of typical all or nothing approach.
Node folding priorities are tracked by a binary heap instead of typical linked list.

The output image is better at preserving textures of the original than Gimp, though it obviously depends on the input image. This particular frog image has the color bar added at the top specifically to throw off some early truncation algorithms, which Gimp is suseptible to. <lang c>typedef struct oct_node_t oct_node_t, *oct_node; struct oct_node_t{ /* sum of all colors represented by this node. 64 bit in case of HUGE image */ uint64_t r, g, b; int count, heap_idx; oct_node kids[8], parent; unsigned char n_kids, kid_idx, flags, depth; };

/* cmp function that decides the ordering in the heap. This is how we determine

  which octree node to fold next, the heart of the algorithm. */

inline int cmp_node(oct_node a, oct_node b) { if (a->n_kids < b->n_kids) return -1; if (a->n_kids > b->n_kids) return 1;

int ac = a->count * (1 + a->kid_idx) >> a->depth; int bc = b->count * (1 + b->kid_idx) >> b->depth; return ac < bc ? -1 : ac > bc; }

/* adding a color triple to octree */ oct_node node_insert(oct_node root, unsigned char *pix) {

define OCT_DEPTH 8

/* 8: number of significant bits used for tree. It's probably good enough for most images to use a value of 5. This affects how many nodes eventually end up in the tree and heap, thus smaller values helps with both speed and memory. */

unsigned char i, bit, depth = 0; for (bit = 1 << 7; ++depth < OCT_DEPTH; bit >>= 1) { i = !!(pix[1] & bit) * 4 + !!(pix[0] & bit) * 2 + !!(pix[2] & bit); if (!root->kids[i]) root->kids[i] = node_new(i, depth, root);

root = root->kids[i]; }

root->r += pix[0]; root->g += pix[1]; root->b += pix[2]; root->count++; return root; }

/* remove a node in octree and add its count and colors to parent node. */ oct_node node_fold(oct_node p) { if (p->n_kids) abort(); oct_node q = p->parent; q->count += p->count;

q->r += p->r; q->g += p->g; q->b += p->b; q->n_kids --; q->kids[p->kid_idx] = 0; return q; }

/* traverse the octree just like construction, but this time we replace the pixel

  color with color stored in the tree node */

void color_replace(oct_node root, unsigned char *pix) { unsigned char i, bit;

for (bit = 1 << 7; bit; bit >>= 1) { i = !!(pix[1] & bit) * 4 + !!(pix[0] & bit) * 2 + !!(pix[2] & bit); if (!root->kids[i]) break; root = root->kids[i]; }

pix[0] = root->r; pix[1] = root->g; pix[2] = root->b; }

/* Building an octree and keep leaf nodes in a bin heap. Afterwards remove first node

  in heap and fold it into its parent node (which may now be added to heap), until heap
  contains required number of colors. */

void color_quant(image im, int n_colors) { int i; unsigned char *pix = im->pix; node_heap heap = { 0, 0, 0 };

oct_node root = node_new(0, 0, 0), got; for (i = 0; i < im->w * im->h; i++, pix += 3) heap_add(&heap, node_insert(root, pix));

while (heap.n > n_colors + 1) heap_add(&heap, node_fold(pop_heap(&heap)));

double c; for (i = 1; i < heap.n; i++) { got = heap.buf[i]; c = got->count; got->r = got->r / c + .5; got->g = got->g / c + .5; got->b = got->b / c + .5; printf("%2d | %3llu %3llu %3llu (%d pixels)\n", i, got->r, got->g, got->b, got->count); }

for (i = 0, pix = im->pix; i < im->w * im->h; i++, pix += 3) color_replace(root, pix);

node_free(); free(heap.buf); }</lang>

OCaml

Here we use a simplified method inspired from this paper: www.leptonica.com/papers/mediancut.pdf

<lang ocaml>let rem_from rem from =

 List.filter ((<>) rem) from

let float_rgb (r,g,b) = (* prevents int overflow *)

 (float r, float g, float b)

let round x =

 int_of_float (floor (x +. 0.5))

let int_rgb (r,g,b) =

 (round r, round g, round b)

let rgb_add (r1,g1,b1) (r2,g2,b2) =

 (r1 +. r2,
  g1 +. g2,
  b1 +. b2)

let rgb_mean px_list =

 let n = float (List.length px_list) in
 let r, g, b = List.fold_left rgb_add (0.0, 0.0, 0.0) px_list in
 (r /. n, g /. n, b /. n)

let extrems lst =

 let min_rgb = (infinity, infinity, infinity)
 and max_rgb = (neg_infinity, neg_infinity, neg_infinity) in
 List.fold_left (fun ((sr,sg,sb), (mr,mg,mb)) (r,g,b) ->
   ((min sr r), (min sg g), (min sb b)),
   ((max mr r), (max mg g), (max mb b))
 ) (min_rgb, max_rgb) lst

let volume_and_dims lst =

 let (sr,sg,sb), (br,bg,bb) = extrems lst in
 let dr, dg, db = (br -. sr), (bg -. sg), (bb -. sb) in
 (dr *. dg *. db),
 (dr, dg, db)

let make_cluster pixel_list =

 let vol, dims = volume_and_dims pixel_list in
 let len = float (List.length pixel_list) in
 (rgb_mean pixel_list, len *. vol, dims, pixel_list)

type axis = R | G | B let largest_axis (r,g,b) =

 match compare r g, compare r b with
 | 1, 1 -> R
 | -1, 1 -> G
 | 1, -1 -> B
 | _ ->
     match compare g b with
     | 1 -> G
     | _ -> B

let subdivide ((mr,mg,mb), n_vol_prod, vol, pixels) =

 let part_func =
   match largest_axis vol with
   | R -> (fun (r,_,_) -> r < mr)
   | G -> (fun (_,g,_) -> g < mg)
   | B -> (fun (_,_,b) -> b < mb)
 in
 let px1, px2 = List.partition part_func pixels in
 (make_cluster px1, make_cluster px2)

let color_quant img n =

 let width, height = get_dims img in
 let clusters =
   let lst = ref [] in
   for x = 0 to pred width do
     for y = 0 to pred height do
       let rgb = float_rgb (get_pixel_unsafe img x y) in
       lst := rgb :: !lst
     done;
   done;
   ref [make_cluster !lst]
 in
 while (List.length !clusters) < n do
   let dumb = (0.0,0.0,0.0) in
   let unused = (dumb, neg_infinity, dumb, []) in
   let select ((_,v1,_,_) as c1) ((_,v2,_,_) as c2) =
     if v1 > v2 then c1 else c2
   in
   let cl = List.fold_left (fun c1 c2 -> select c1 c2) unused !clusters in
   let cl1, cl2 = subdivide cl in
   clusters := cl1 :: cl2 :: (rem_from cl !clusters)
 done;
 let module PxMap = Map.Make
   (struct type t = float * float * float let compare = compare end) in
 let m =
   List.fold_left (fun m (mean, _, _, pixel_list) ->
     let int_mean = int_rgb mean in
     List.fold_left (fun m px -> PxMap.add px int_mean m) m pixel_list
   ) PxMap.empty !clusters
 in
 let res = new_img ~width ~height in
 for y = 0 to pred height do
   for x = 0 to pred width do
     let rgb = float_rgb (get_pixel_unsafe img x y) in
     let mean_rgb = PxMap.find rgb m in
     put_pixel_unsafe res mean_rgb x y;
   done;
 done;
 (res)</lang>

J

Here, we use a simplistic averaging technique to build an initial set of colors and then use k-means clustering to refine them.

<lang j>kmcL=:4 :0

 C=. /:~ 256 #.inv ,y  NB. colors
 G=. x (i.@] <.@* %) #C  NB. groups (initial)
 Q=. _  NB. quantized list of colors (initial
 whilst.-. Q-:&<.&(x&*)Q0 do.
   Q0=. Q
   Q=. /:~C (+/ % #)/.~ G
   G=. (i. <./)"1 C +/&.:*: .- |:Q
 end.Q

)</lang>

The left argument is the number of colors desired.

The right argument is the image, with pixels represented as bmp color integers (base 256 numbers).

The result is the colors represented as pixel triples (blue, green, red). They are shown here as fractional numbers, but they should be either rounded to the nearest integer in the range 0..255 (and possibly converted back to bmp integer form) or scaled so they are floating point triples in the range 0..1.