| 1 |
/* |
| 2 |
|
| 3 |
#ifndef lint |
| 4 |
static char SCCSid[] = "$SunId$ LBL"; |
| 5 |
#endif |
| 6 |
CUT - Median bisection (k-d tree) algorithm for color image quantization |
| 7 |
|
| 8 |
flaw: doesn't always create as many colors as requested |
| 9 |
because the binary subdivision is always done in a balanced tree and when a box |
| 10 |
contains only one color it's unsplittable |
| 11 |
|
| 12 |
Paul Heckbert |
| 13 |
*/ |
| 14 |
|
| 15 |
#include "ciq.h" |
| 16 |
|
| 17 |
struct index {short *o; int freq;}; |
| 18 |
/* offset and pixel count for a box */ |
| 19 |
|
| 20 |
struct pail {struct plum *first; int freak;}; |
| 21 |
/* each pail is a bucket containing plums */ |
| 22 |
|
| 23 |
struct plum {struct plum *next; int code;}; |
| 24 |
/* a color in a pail's linked list */ |
| 25 |
|
| 26 |
static short off[len]; |
| 27 |
/* color codes: offsets of colors in hist array */ |
| 28 |
|
| 29 |
|
| 30 |
makecm(nw,na) /* subdivide colorspace to generate a colormap*/ |
| 31 |
int nw,*na; /* number of colors wanted, actual */ |
| 32 |
{ |
| 33 |
int i,m,n,freq,weight,sr,sg,sb; |
| 34 |
short *o; |
| 35 |
struct index index[257]; /* pointers into off */ |
| 36 |
/* index[i].o is pointer to code for first color in box i */ |
| 37 |
|
| 38 |
for (o=off, *na=i=0; i<len; i++) |
| 39 |
if (hist[i]) { |
| 40 |
(*na)++; |
| 41 |
*o++ = i; /* initialize offsets */ |
| 42 |
} |
| 43 |
index[0].o = off; /* start with one big box which contains all *na colors */ |
| 44 |
index[0].freq = xmax*ymax; /* and (sum of hist[i]) pixels */ |
| 45 |
index[1].o = off+*na; |
| 46 |
|
| 47 |
/* breadth-first binary subdivision: try to make nw boxes */ |
| 48 |
for (m=1; m<nw; m<<=1) { |
| 49 |
/*printf("%d ",m);*/ |
| 50 |
index[m<<1].o = index[m].o; |
| 51 |
for (i=m-1; i>=0; --i) splitbox(index+i,index+(i<<1)); |
| 52 |
} |
| 53 |
for (n=i=0; i<m && n<nw; i++) if (index[i+1].o>index[i].o) { |
| 54 |
sr = sg = sb = (freq = index[i].freq)>>1; |
| 55 |
for (o=index[i].o; o<index[i+1].o; o++) { |
| 56 |
weight = hist[*o]; |
| 57 |
/* find weighted average of colors in box i */ |
| 58 |
sr += red(*o)*weight; |
| 59 |
sg += gre(*o)*weight; |
| 60 |
sb += blu(*o)*weight; |
| 61 |
} |
| 62 |
color[0][n] = sr/freq; |
| 63 |
color[1][n] = sg/freq; |
| 64 |
color[2][n] = sb/freq; |
| 65 |
/* printf("color[%d]=(%3d,%3d,%3d) freq=%d\n", |
| 66 |
n,color[0][n],color[1][n],color[2][n],freq); */ |
| 67 |
n++; |
| 68 |
} |
| 69 |
/*printf("[%d empty boxes] from %d to %d colors\n",nw-n,*na,n);*/ |
| 70 |
return n; |
| 71 |
} |
| 72 |
|
| 73 |
splitbox(ii,io) /* split a box in two: half of the pixels from */ |
| 74 |
struct index *ii,*io; /* box ii will go into io, the other half into io+1 */ |
| 75 |
{ |
| 76 |
register short *o,*o1,*o2; |
| 77 |
register int shift,count,k,freq,r,g,b,r1,g1,b1,r2,g2,b2; |
| 78 |
struct pail pail[32],*p,*q; /* buckets for sorting colors */ |
| 79 |
struct plum plum[len],*pl=plum,*u; /* colors */ |
| 80 |
|
| 81 |
o1 = ii[0].o; |
| 82 |
o2 = ii[1].o; |
| 83 |
freq = ii[0].freq; /* number of pixels with color in this box */ |
| 84 |
r1 = g1 = b1 = 256; r2 = g2 = b2 = -1; |
| 85 |
for (o=o1; o<o2; o++) { |
| 86 |
r = red(*o); |
| 87 |
g = gre(*o); |
| 88 |
b = blu(*o); |
| 89 |
if (r<r1) r1 = r; if (r>r2) r2 = r; /* find min&max r, g, and b */ |
| 90 |
if (g<g1) g1 = g; if (g>g2) g2 = g; |
| 91 |
if (b<b1) b1 = b; if (b>b2) b2 = b; |
| 92 |
} |
| 93 |
|
| 94 |
/* adaptive partitioning: find dominant dimension */ |
| 95 |
shift = r2-r1>g2-g1?r2-r1>=b2-b1?10:0:g2-g1>=b2-b1?5:0; |
| 96 |
/* printf("splitbox %3d-%3d %3dx%3dx%3d %c ", |
| 97 |
o1-off,o2-off,r2-r1,g2-g1,b2-b1,shift==10?'R':shift==5?'G':'B'); */ |
| 98 |
|
| 99 |
for (p=pail, k=0; k<32; k++, p++) { /* start with empty pails */ |
| 100 |
p->first = 0; |
| 101 |
p->freak = 0; |
| 102 |
} |
| 103 |
for (o=o1; o<o2; o++) { |
| 104 |
p = pail+(*o>>shift&31); /* find appropriate pail */ |
| 105 |
pl->code = *o; |
| 106 |
pl->next = p->first; |
| 107 |
p->first = pl++; /* add new plum to pail */ |
| 108 |
p->freak += hist[*o]; |
| 109 |
} |
| 110 |
|
| 111 |
/* find median along dominant dimension */ |
| 112 |
for (count=freq>>1, p=pail; count>0; p++) count -= p->freak; |
| 113 |
|
| 114 |
if (p>pail && count+(p-1)->freak<-count) /* back up one? */ |
| 115 |
count += (--p)->freak; |
| 116 |
|
| 117 |
io[1].o = o1; /* in case of empty box */ |
| 118 |
for (o=o1, q=pail; o<o2;) { /* dump the pails to reorder colors in box */ |
| 119 |
for (u=q->first; u; u=u->next) *o++ = u->code; |
| 120 |
if (++q==p) io[1].o = o; /* place partition at offset o */ |
| 121 |
} |
| 122 |
io[0].o = o1; |
| 123 |
io[0].freq = (freq>>1)-count; |
| 124 |
io[1].freq = (freq+1>>1)+count; |
| 125 |
/* printf(" at %3d %d=%d+%d\n",io[1].o-off,freq,io[0].freq,io[1].freq); */ |
| 126 |
} |
