Color quantization: Difference between revisions
Line 522: | Line 522: | ||
[[file:Compare_16_Quantum_frog_histograms_PureBasic.png|histogram (external application)|thumb|200px]] |
[[file:Compare_16_Quantum_frog_histograms_PureBasic.png|histogram (external application)|thumb|200px]] |
||
<lang PureBasic> |
<lang PureBasic> |
||
; ColorQuantization.pb |
|||
Structure bestA_ ; table for our histogram |
|||
nn.i ; 16,32,... |
|||
rc.i ; red count within (0,1,...,255)/(number of colors) |
|||
gc.i ; green count within (0,1,...,255)/(number of colors) |
|||
bc.i ; blue count within (0,1,...,255)/(number of colors) |
|||
EndStructure |
|||
; these two functions appear to be rather self-explanatory |
|||
</lang> |
|||
UsePNGImageDecoder() |
|||
UsePNGImageEncoder() |
|||
Procedure.i ColorQuantization(Filename$,ncol) |
|||
'''Sample Output''' |
|||
Protected x,y,c |
|||
<pre> |
|||
; load our original image or leave the procedure |
|||
</pre> |
|||
If not LoadImage(0,Filename$) :ProcedureReturn 0:endif |
|||
; we are not going to actually draw on the original image... |
|||
; but we need to use the drawing library to load up |
|||
; the pixel information into our arrays... |
|||
; if we can't do that, what's the point of going any further? |
|||
; so then we would be wise to just leave the procedure [happy fred?] |
|||
If not StartDrawing(ImageOutput(0)):ProcedureReturn 0:endif |
|||
iw=ImageWidth(0) |
|||
ih=ImageHeight(0) |
|||
dim cA(iw,ih) ; color array to hold at a given (x,y) |
|||
dim rA(iw,ih) ; red array to hold at a given (x,y) |
|||
dim gA(iw,ih) ; green array to hold at a given (x,y) |
|||
dim bA(iw,ih) ; blue array to hold at a given (x,y) |
|||
dim tA(iw,ih) ; temp array to hold at a given (x,y) |
|||
; map each pixel from the original image to our arrays |
|||
; don't overrun the ranges ie. use {ih-1,iw-1} |
|||
for y=0 to ih-1 |
|||
for x=0 to iw-1 |
|||
c = Point(x,y) |
|||
cA(x,y)=c |
|||
rA(x,y)=Red(c) |
|||
gA(x,y)=Green(c) |
|||
bA(x,y)=Blue(c) |
|||
next |
|||
next |
|||
StopDrawing() ; don't forget to... StopDrawing() |
|||
N=ih*iw |
|||
; N is the total number if pixels |
|||
if not N:ProcedureReturn 0:endif ; to avoid a division by zero |
|||
; stuctured array ie. a table to hold the frequency distribution |
|||
dim bestA.bestA_(ncol) |
|||
; the "best" red,green,blue based upon frequency |
|||
dim rbestA(ncol/3) |
|||
dim gbestA(ncol/3) |
|||
dim bbestA(ncol/3) |
|||
; split the (0..255) range up |
|||
xoff=256/ncol ;256/16=16 |
|||
xrng=xoff ;xrng=16 |
|||
; store these values in our table: bestA(i)\nn= 16,32,... |
|||
for i=1 to ncol |
|||
xrng+xoff |
|||
bestA(i)\nn=xrng |
|||
next |
|||
; scan by row [y] |
|||
for y=0 to ih-1 |
|||
; scan by col [x] |
|||
for x=0 to iw-1 |
|||
; retrieve the rgb values from each pixel |
|||
r=rA(x,y) |
|||
g=gA(x,y) |
|||
b=bA(x,y) |
|||
; sum up the numbers that fall within our subdivisions of (0..255) |
|||
for i=1 to ncol |
|||
if r>=bestA(i)\nn and r<bestA(i+1)\nn:bestA(i)\rc+1:endif |
|||
if g>=bestA(i)\nn and g<bestA(i+1)\nn:bestA(i)\gc+1:endif |
|||
if b>=bestA(i)\nn and b<bestA(i+1)\nn:bestA(i)\bc+1:endif |
|||
next |
|||
next |
|||
next |
|||
; option and type to: Sort our Structured Array |
|||
opt=#PB_Sort_Descending |
|||
typ=#PB_Sort_Integer |
|||
; sort to get most frequent reds |
|||
off=OffsetOf(bestA_\rc) |
|||
SortStructuredArray(bestA(),opt, off, typ,1, ncol) |
|||
; save the best [ for number of colors =16 this is int(16/3)=5 ] reds |
|||
for i=1 to ncol/3 |
|||
rbestA(i)=bestA(i)\nn |
|||
next |
|||
; sort to get most frequent greens |
|||
off=OffsetOf(bestA_\gc) |
|||
SortStructuredArray(bestA(),opt, off, typ,1, ncol) |
|||
; save the best [ for number of colors =16 this is int(16/3)=5 ] greens |
|||
for i=1 to ncol/3 |
|||
gbestA(i)=bestA(i)\nn |
|||
next |
|||
; sort to get most frequent blues |
|||
off=OffsetOf(bestA_\bc) |
|||
SortStructuredArray(bestA(),opt, off, typ,1, ncol) |
|||
; save the best [ for number of colors =16 this is int(16/3)=5 ] blues |
|||
for i=1 to ncol/3 |
|||
bbestA(i)=bestA(i)\nn |
|||
next |
|||
; reset the best low value to 15 and high value to 240 |
|||
; this helps to ensure there is some contrast when the statistics bunch up |
|||
; ie. when a single color tends to predominate... such as perhaps green? |
|||
rbestA(1)=15:rbestA(ncol/3)=240 |
|||
gbestA(1)=15:gbestA(ncol/3)=240 |
|||
bbestA(1)=15:bbestA(ncol/3)=240 |
|||
; make a copy of our original image or leave the procedure |
|||
If not CopyImage(0,1) :ProcedureReturn 0:endif |
|||
; draw on that copy of our original image or leave the procedure |
|||
If not StartDrawing(ImageOutput(1)):ProcedureReturn 0:endif |
|||
for y=0 to ih-1 |
|||
for x=0 to iw-1 |
|||
c = Point(x,y) |
|||
; get the rgb value from our arrays |
|||
rt=rA(x,y) |
|||
gt=gA(x,y) |
|||
bt=bA(x,y) |
|||
; given a particular red value say 123 at point x,y |
|||
; which of our rbestA(i's) is closest? |
|||
; then for green and blue? |
|||
; ============================== |
|||
r=255 |
|||
for i=1 to ncol/3 |
|||
rdiff=abs(rbestA(i)-rt) |
|||
if rdiff<=r:ri=i:r=rdiff:endif |
|||
next |
|||
g=255 |
|||
for i=1 to ncol/3 |
|||
gdiff=abs(gbestA(i)-gt) |
|||
if gdiff<=g:gi=i:g=gdiff:endif |
|||
next |
|||
b=255 |
|||
for i=1 to ncol/3 |
|||
bdiff=abs(bbestA(i)-bt) |
|||
if bdiff<=b:bi=i:b=bdiff:endif |
|||
next |
|||
; ============================== |
|||
; get the color value so we can plot it at that pixel |
|||
Color=RGB(rbestA(ri),gbestA(gi),bbestA(bi)) |
|||
; plot it at that pixel |
|||
Plot(x,y,Color) |
|||
; save that info to tA(x,y) for our comparison image |
|||
tA(x,y)=Color |
|||
next |
|||
next |
|||
StopDrawing() ; don't forget to... StopDrawing() |
|||
; create a comparison image of our original vs 16-color or leave the procedure |
|||
If not CreateImage(2,iw*2,ih) :ProcedureReturn 0:endif |
|||
; draw on that image both our original image and our 16-color image or leave the procedure |
|||
If not StartDrawing(ImageOutput(2)):ProcedureReturn 0:endif |
|||
; plot original image |
|||
; 0,0 .... 511,0 |
|||
; . |
|||
; . |
|||
; 511,0 .. 511,511 |
|||
for y=0 to ih-1 |
|||
for x=0 to iw-1 |
|||
c = cA(x,y) |
|||
Plot(x,y,c) |
|||
next |
|||
next |
|||
; plot 16-color image to the right of original image |
|||
; 512,0 .... 1023,0 |
|||
; . |
|||
; . |
|||
; 512,511 .. 1023,511 |
|||
for y=0 to ih-1 |
|||
for x=0 to iw-1 |
|||
c = tA(x,y) |
|||
Plot(x+iw,y,c) |
|||
next |
|||
next |
|||
StopDrawing() ; don't forget to... StopDrawing() |
|||
; save the single 16-color image |
|||
SaveImage(1, "_single_"+str(ncol)+"_"+Filename$,#PB_ImagePlugin_PNG ) |
|||
; save the comparison image |
|||
SaveImage(2, "_compare_"+str(ncol)+"_"+Filename$,#PB_ImagePlugin_PNG ) |
|||
ProcedureReturn 1 |
|||
EndProcedure |
|||
ColorQuantization("Quantum_frog.png",16) |
|||
</lang> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
Revision as of 23:27, 18 August 2012
You are encouraged to solve this task according to the task description, using any language you may know.
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>
D
This code retains some of the style of the original OCaml code. <lang d>import core.stdc.stdio, std.stdio, std.ascii, std.algorithm,
std.typecons, std.math, std.range, std.conv;
final class Image {
int w, h; ubyte[] pix;
void allocate(in int nr, in int nc) pure nothrow { w = nc; h = nr; this.pix.length = 3 * this.w * this.h; }
void loadPPM6(in string fileName) { static int read_num(FILE* f) nothrow { int n; while (!fscanf(f, "%d ", &n)) { if ((n = fgetc(f)) == '#') { while ((n = fgetc(f)) != '\n') if (n == EOF) return 0; } else return 0; } return n; }
auto fin = fopen((fileName ~ '\0').ptr, "rb"); if (!fin) goto bail;
if (fgetc(fin) != 'P' || fgetc(fin) != '6' || !isWhite(fgetc(fin))) goto bail;
immutable int nc = read_num(fin); immutable int nr = read_num(fin); immutable int maxval = read_num(fin); if (nc <= 0 || nr <= 0 || maxval <= 0) goto bail; allocate(nr, nc);
immutable count = fread(this.pix.ptr, 1, 3 * nc * nr, fin); if (count != 3 * nc * nr) throw new Exception("Wrong number of items read.");
bail: // scope(exit) is better if (fin) fclose(fin); }
void savePPM6(in string fileName) const in { assert(!fileName.empty); assert(this.w > 0 && this.h > 0 && pix.length == (3 * this.w * this.h), "Not correct image."); } body { auto fout = fopen((fileName ~ '\0').ptr, "wb"); if (fout == null) throw new Exception("File can't be opened."); fprintf(fout, "P6\n%d %d\n255\n", this.w, this.h); immutable count = fwrite(this.pix.ptr, 1, 3 * this.w * this.h, fout); if (count != 3 * this.w * this.h) new Exception("Wrong number of items written."); fclose(fout); }
}
alias Tuple!(float,"r", float,"g", float,"b") Col; alias Tuple!(Col, float, Col, Col[]) Cluster; enum Axis { R, G, B }
// ubyte[3] causes slow heap allocations alias Tuple!(ubyte, ubyte, ubyte) Ubyte3;
int round(in float x) /*pure*/ nothrow {
return cast(int)floor(x + 0.5); // not pure
}
Ubyte3 roundUbyteRGB(in Col c) /*pure*/ nothrow {
return tuple(cast(ubyte)round(c.r), cast(ubyte)round(c.g), cast(ubyte)round(c.b));
}
Col meanRGB(Col[] pxList) pure nothrow {
static Col addRGB(in Col c1, in Col c2) pure nothrow { return Col(c1.r + c2.r, c1.g + c2.g, c1.b + c2.b); } immutable Col tot = reduce!addRGB(Col(0,0,0), pxList); immutable int n = walkLength(pxList); return Col(tot.r / n, tot.g / n, tot.b / n);
}
Tuple!(Col, Col) extrems(Col[] lst) pure nothrow {
immutable minRGB = Col(float.infinity, float.infinity, float.infinity); immutable maxRGB = Col(-float.infinity, -float.infinity, -float.infinity); static Col f1(in Col c1, in Col c2) pure nothrow { return Col(min(c1.r, c2.r), min(c1.g, c2.g), min(c1.b, c2.b)); } static Col f2(in Col c1, in Col c2) pure nothrow { return Col(max(c1.r, c2.r), max(c1.g, c2.g), max(c1.b, c2.b)); } return typeof(return)(reduce!f1(minRGB, lst), reduce!f2(maxRGB, lst));
}
Tuple!(float, Col) volumeAndDims(Col[] lst) pure nothrow {
immutable e = extrems(lst); immutable Col r = Col(e[1].r - e[0].r, e[1].g - e[0].g, e[1].b - e[0].b); return typeof(return)(r.r * r.g * r.b, r);
}
Cluster makeCluster(Col[] pixel_list) pure nothrow {
immutable vol_dims = volumeAndDims(pixel_list); immutable int len = pixel_list.length; return Cluster(meanRGB(pixel_list), len * vol_dims[0], vol_dims[1], pixel_list);
}
Axis largestAxis(in Col c) pure nothrow {
static int fcmp(in float a, in float b) pure nothrow { return (a > b) ? 1 : (a < b ? -1 : 0); } immutable int r1 = fcmp(c.r, c.g); immutable int r2 = fcmp(c.r, c.b); if (r1 == 1 && r2 == 1) return Axis.R; if (r1 == -1 && r2 == 1) return Axis.G; if (r1 == 1 && r2 == -1) return Axis.B; return (fcmp(c.g, c.b) == 1) ? Axis.G : Axis.B;
}
Tuple!(Cluster, Cluster) subdivide(in Col c, in float nVolProd,
in Col vol, Col[] pixels)
pure /*nothrow*/ {
bool delegate(Col c) partFunc; final switch (largestAxis(vol)) { case Axis.R: partFunc = c1 => c1.r < c.r; break; case Axis.G: partFunc = c1 => c1.g < c.g; break; case Axis.B: partFunc = c1 => c1.b < c.b; break; } Col[] px2 = partition!partFunc(pixels); Col[] px1 = pixels[0 .. $ - px2.length]; return typeof(return)(makeCluster(px1), makeCluster(px2));
}
Image colorQuantize(in Image img, in int n) /*pure nothrow*/ {
immutable int width = img.w; immutable int height = img.h;
auto cols = new Col[width * height]; foreach (i, ref c; cols) c = Col(img.pix[i * 3 + 0], img.pix[i * 3 + 1], img.pix[i * 3 + 2]); Cluster[] clusters = [makeCluster(cols)];
immutable Col dumb = Col(0.0, 0.0, 0.0); Cluster unused = Cluster(dumb, -float.infinity, dumb, (Col[]).init);
while (clusters.length < n) { Cluster cl = reduce!((c1,c2) => c1[1] > c2[1] ? c1 : c2) (unused, clusters); clusters = [subdivide(cl.tupleof).tupleof] ~ remove!(c => c == cl, SwapStrategy.unstable)(clusters); }
static uint ubyte3ToUint(in Ubyte3 c) pure nothrow { uint r; r |= c[0]; r |= c[1] << 8; r |= c[2] << 16; return r; }
uint[uint] pixMap; // faster than Ubyte3[Ubyte3] ubyte[4] u4a, u4b; foreach (cluster; clusters) { immutable ubyteMean = ubyte3ToUint(roundUbyteRGB(cluster[0])); foreach (col; cluster[3]) pixMap[ubyte3ToUint(roundUbyteRGB(col))] = ubyteMean; }
auto result = new Image; result.allocate(height, width);
static Ubyte3 uintToUbyte3(in uint c) pure nothrow { Ubyte3 r = void; r[0] = c & 0b1111_1111; r[1] = (c >> 8) & 0b1111_1111; r[2] = (c >> 16) & 0b1111_1111; return r; }
foreach (i; 0 .. height * width) { immutable Ubyte3 u3a = Ubyte3(img.pix[i * 3 + 0], img.pix[i * 3 + 1], img.pix[i * 3 + 2]); immutable Ubyte3 c = uintToUbyte3(pixMap[ubyte3ToUint(u3a)]); result.pix[i * 3 + 0] = c[0]; result.pix[i * 3 + 1] = c[1]; result.pix[i * 3 + 2] = c[2]; }
return result;
}
void main(in string[] args) {
string fileName; int nCols; if (args.length == 3) { fileName = args[1]; nCols = to!int(args[2]); } else if (args.length == 1) { fileName = "quantum_frog.ppm"; nCols = 16; } else { writeln("Usage: color_quantization image.ppm ncolors"); return; }
auto im = new Image; im.loadPPM6(fileName); const imq = colorQuantize(im, nCols); imq.savePPM6("out.ppm");
}</lang>
Mathematica
<lang Mathematica>ColorQuantize[Import["http://rosettacode.org/mw/images/3/3f/Quantum_frog.png"],16,Dithering->False]</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.
<lang j> 16 kmcL img 7.52532 22.3347 0.650468 8.20129 54.4678 0.0326828 33.1132 69.8148 0.622265 54.2232 125.682 2.67713 56.7064 99.5008 3.04013 61.2135 136.42 4.2015 68.1246 140.576 6.37512 74.6006 143.606 7.57854 78.9101 150.792 10.2563 89.5873 148.621 14.6202 98.9523 154.005 25.7583 114.957 159.697 47.6423 145.816 178.136 33.8845 164.969 199.742 67.0467 179.849 207.594 109.973 209.229 221.18 204.513</lang>
PureBasic
<lang PureBasic>
- ColorQuantization.pb
Structure bestA_ ; table for our histogram nn.i ; 16,32,... rc.i ; red count within (0,1,...,255)/(number of colors) gc.i ; green count within (0,1,...,255)/(number of colors) bc.i ; blue count within (0,1,...,255)/(number of colors) EndStructure
- these two functions appear to be rather self-explanatory
UsePNGImageDecoder() UsePNGImageEncoder()
Procedure.i ColorQuantization(Filename$,ncol) Protected x,y,c
- load our original image or leave the procedure
If not LoadImage(0,Filename$) :ProcedureReturn 0:endif
- we are not going to actually draw on the original image...
- but we need to use the drawing library to load up
- the pixel information into our arrays...
- if we can't do that, what's the point of going any further?
- so then we would be wise to just leave the procedure [happy fred?]
If not StartDrawing(ImageOutput(0)):ProcedureReturn 0:endif
iw=ImageWidth(0) ih=ImageHeight(0)
dim cA(iw,ih) ; color array to hold at a given (x,y) dim rA(iw,ih) ; red array to hold at a given (x,y) dim gA(iw,ih) ; green array to hold at a given (x,y) dim bA(iw,ih) ; blue array to hold at a given (x,y) dim tA(iw,ih) ; temp array to hold at a given (x,y)
- map each pixel from the original image to our arrays
- don't overrun the ranges ie. use {ih-1,iw-1}
for y=0 to ih-1
for x=0 to iw-1 c = Point(x,y) cA(x,y)=c rA(x,y)=Red(c) gA(x,y)=Green(c) bA(x,y)=Blue(c) next
next
StopDrawing() ; don't forget to... StopDrawing()
N=ih*iw
- N is the total number if pixels
if not N:ProcedureReturn 0:endif ; to avoid a division by zero
- stuctured array ie. a table to hold the frequency distribution
dim bestA.bestA_(ncol)
- the "best" red,green,blue based upon frequency
dim rbestA(ncol/3) dim gbestA(ncol/3) dim bbestA(ncol/3)
- split the (0..255) range up
xoff=256/ncol ;256/16=16 xrng=xoff ;xrng=16
- store these values in our table
- bestA(i)\nn= 16,32,...
for i=1 to ncol xrng+xoff bestA(i)\nn=xrng next
- scan by row [y]
for y=0 to ih-1
- scan by col [x]
for x=0 to iw-1
- retrieve the rgb values from each pixel
r=rA(x,y) g=gA(x,y) b=bA(x,y)
- sum up the numbers that fall within our subdivisions of (0..255)
for i=1 to ncol if r>=bestA(i)\nn and r<bestA(i+1)\nn:bestA(i)\rc+1:endif if g>=bestA(i)\nn and g<bestA(i+1)\nn:bestA(i)\gc+1:endif if b>=bestA(i)\nn and b<bestA(i+1)\nn:bestA(i)\bc+1:endif next next next
- option and type to
- Sort our Structured Array
opt=#PB_Sort_Descending typ=#PB_Sort_Integer
- sort to get most frequent reds
off=OffsetOf(bestA_\rc) SortStructuredArray(bestA(),opt, off, typ,1, ncol)
- save the best [ for number of colors =16 this is int(16/3)=5 ] reds
for i=1 to ncol/3 rbestA(i)=bestA(i)\nn next
- sort to get most frequent greens
off=OffsetOf(bestA_\gc) SortStructuredArray(bestA(),opt, off, typ,1, ncol)
- save the best [ for number of colors =16 this is int(16/3)=5 ] greens
for i=1 to ncol/3 gbestA(i)=bestA(i)\nn next
- sort to get most frequent blues
off=OffsetOf(bestA_\bc) SortStructuredArray(bestA(),opt, off, typ,1, ncol)
- save the best [ for number of colors =16 this is int(16/3)=5 ] blues
for i=1 to ncol/3 bbestA(i)=bestA(i)\nn next
- reset the best low value to 15 and high value to 240
- this helps to ensure there is some contrast when the statistics bunch up
- ie. when a single color tends to predominate... such as perhaps green?
rbestA(1)=15:rbestA(ncol/3)=240 gbestA(1)=15:gbestA(ncol/3)=240 bbestA(1)=15:bbestA(ncol/3)=240
- make a copy of our original image or leave the procedure
If not CopyImage(0,1) :ProcedureReturn 0:endif
- draw on that copy of our original image or leave the procedure
If not StartDrawing(ImageOutput(1)):ProcedureReturn 0:endif
for y=0 to ih-1 for x=0 to iw-1 c = Point(x,y)
- get the rgb value from our arrays
rt=rA(x,y) gt=gA(x,y) bt=bA(x,y)
- given a particular red value say 123 at point x,y
- which of our rbestA(i's) is closest?
- then for green and blue?
- ==============================
r=255 for i=1 to ncol/3 rdiff=abs(rbestA(i)-rt) if rdiff<=r:ri=i:r=rdiff:endif next
g=255 for i=1 to ncol/3 gdiff=abs(gbestA(i)-gt) if gdiff<=g:gi=i:g=gdiff:endif next
b=255 for i=1 to ncol/3 bdiff=abs(bbestA(i)-bt) if bdiff<=b:bi=i:b=bdiff:endif next
- ==============================
- get the color value so we can plot it at that pixel
Color=RGB(rbestA(ri),gbestA(gi),bbestA(bi))
- plot it at that pixel
Plot(x,y,Color)
- save that info to tA(x,y) for our comparison image
tA(x,y)=Color
next next StopDrawing() ; don't forget to... StopDrawing()
- create a comparison image of our original vs 16-color or leave the procedure
If not CreateImage(2,iw*2,ih) :ProcedureReturn 0:endif
- draw on that image both our original image and our 16-color image or leave the procedure
If not StartDrawing(ImageOutput(2)):ProcedureReturn 0:endif
- plot original image
- 0,0 .... 511,0
- .
- .
- 511,0 .. 511,511
for y=0 to ih-1
for x=0 to iw-1 c = cA(x,y) Plot(x,y,c) next next
- plot 16-color image to the right of original image
- 512,0 .... 1023,0
- .
- .
- 512,511 .. 1023,511
for y=0 to ih-1
for x=0 to iw-1 c = tA(x,y) Plot(x+iw,y,c) next next
StopDrawing() ; don't forget to... StopDrawing()
- save the single 16-color image
SaveImage(1, "_single_"+str(ncol)+"_"+Filename$,#PB_ImagePlugin_PNG )
- save the comparison image
SaveImage(2, "_compare_"+str(ncol)+"_"+Filename$,#PB_ImagePlugin_PNG ) ProcedureReturn 1 EndProcedure
ColorQuantization("Quantum_frog.png",16)
</lang>
Tcl
<lang tcl>package require Tcl 8.6 package require Tk
proc makeCluster {pixels} {
set rmin [set rmax [lindex $pixels 0 0]] set gmin [set gmax [lindex $pixels 0 1]] set bmin [set bmax [lindex $pixels 0 2]] set rsum [set gsum [set bsum 0]] foreach p $pixels {
lassign $p r g b if {$r<$rmin} {set rmin $r} elseif {$r>$rmax} {set rmax $r} if {$g<$gmin} {set gmin $g} elseif {$g>$gmax} {set gmax $g} if {$b<$bmin} {set bmin $b} elseif {$b>$bmax} {set bmax $b} incr rsum $r incr gsum $g incr bsum $b
} set n [llength $pixels] list [expr {double($n)*($rmax-$rmin)*($gmax-$gmin)*($bmax-$bmin)}] \
[list [expr {$rsum/$n}] [expr {$gsum/$n}] [expr {$bsum/$n}]] \ [list [expr {$rmax-$rmin}] [expr {$gmax-$gmin}] [expr {$bmax-$bmin}]] \ $pixels }
proc colorQuant {img n} {
set width [image width $img] set height [image height $img] # Extract the pixels from the image for {set x 0} {$x < $width} {incr x} {
for {set y 0} {$y < $height} {incr y} { lappend pixels [$img get $x $y] }
} # Divide pixels into clusters for {set cs [list [makeCluster $pixels]]} {[llength $cs] < $n} {} {
set cs [lsort -real -index 0 $cs] lassign [lindex $cs end] score centroid volume pixels lassign $centroid cr cg cb lassign $volume vr vg vb while 1 { set p1 [set p2 {}] if {$vr>$vg && $vr>$vb} { foreach p $pixels { if {[lindex $p 0]<$cr} {lappend p1 $p} {lappend p2 $p} } } elseif {$vg>$vb} { foreach p $pixels { if {[lindex $p 1]<$cg} {lappend p1 $p} {lappend p2 $p} } } else { foreach p $pixels { if {[lindex $p 2]<$cb} {lappend p1 $p} {lappend p2 $p} } } if {[llength $p1] && [llength $p2]} break # Partition failed! Perturb partition point away from the centroid and try again set cr [expr {$cr + 20*rand() - 10}] set cg [expr {$cg + 20*rand() - 10}] set cb [expr {$cb + 20*rand() - 10}] } set cs [lreplace $cs end end [makeCluster $p1] [makeCluster $p2]]
} # Produce map from pixel values to quantized values foreach c $cs {
set centroid [format "#%02x%02x%02x" {*}[lindex $c 1]] foreach p [lindex $c end] { set map($p) $centroid }
} # Remap the source image set newimg [image create photo -width $width -height $height] for {set x 0} {$x < $width} {incr x} {
for {set y 0} {$y < $height} {incr y} { $newimg put $map([$img get $x $y]) -to $x $y }
} return $newimg
}</lang> Demonstration code: <lang tcl>set src [image create photo -file quantum_frog.png] set dst [colorQuant $src 16]
- Save as GIF now that quantization is done, then exit explicitly (no GUI desired)
$dst write quantum_frog_compressed.gif exit</lang>