1 ///////////////////////////////////////////////////////////////////////
2 // C Implementation of Wu's Color Quantizer (v. 2)
3 // (see Graphics Gems vol. II, pp. 126-133)
4 //
5 // Author: Xiaolin Wu
6 // Dept. of Computer Science
7 // Univ. of Western Ontario
8 // London, Ontario N6A 5B7
9 // wu@csd.uwo.ca
10 //
11 // Algorithm: Greedy orthogonal bipartition of RGB space for variance
12 // minimization aided by inclusion-exclusion tricks.
13 // For speed no nearest neighbor search is done. Slightly
14 // better performance can be expected by more sophisticated
15 // but more expensive versions.
16 //
17 // The author thanks Tom Lane at Tom_Lane@G.GP.CS.CMU.EDU for much of
18 // additional documentation and a cure to a previous bug.
19 //
20 // Free to distribute, comments and suggestions are appreciated.
21 ///////////////////////////////////////////////////////////////////////
22
23 ///////////////////////////////////////////////////////////////////////
24 // History
25 // -------
26 // July 2000: C++ Implementation of Wu's Color Quantizer
27 // and adaptation for the FreeImage 2 Library
28 // Author: Herv� Drolon (drolon@infonie.fr)
29 // March 2004: Adaptation for the FreeImage 3 library (port to big endian processors)
30 // Author: Herv� Drolon (drolon@infonie.fr)
31 ///////////////////////////////////////////////////////////////////////
32
33 #include "Quantizers.h"
34 #include "FreeImage.h"
35 #include "Utilities.h"
36
37 ///////////////////////////////////////////////////////////////////////
38
39 // Size of a 3D array : 33 x 33 x 33
40 #define SIZE_3D 35937
41
42 // 3D array indexation
43 #define INDEX(r, g, b) ((r << 10) + (r << 6) + r + (g << 5) + g + b)
44
45 #define MAXCOLOR 256
46
47 // Constructor / Destructor
48
WuQuantizer(FIBITMAP * dib)49 WuQuantizer::WuQuantizer(FIBITMAP *dib) {
50 width = FreeImage_GetWidth(dib);
51 height = FreeImage_GetHeight(dib);
52 pitch = FreeImage_GetPitch(dib);
53 m_dib = dib;
54
55 gm2 = NULL;
56 wt = mr = mg = mb = NULL;
57 Qadd = NULL;
58
59 // Allocate 3D arrays
60 gm2 = (float*)malloc(SIZE_3D * sizeof(float));
61 wt = (LONG*)malloc(SIZE_3D * sizeof(LONG));
62 mr = (LONG*)malloc(SIZE_3D * sizeof(LONG));
63 mg = (LONG*)malloc(SIZE_3D * sizeof(LONG));
64 mb = (LONG*)malloc(SIZE_3D * sizeof(LONG));
65
66 // Allocate Qadd
67 Qadd = (WORD *)malloc(sizeof(WORD) * width * height);
68
69 if(!gm2 || !wt || !mr || !mg || !mb || !Qadd) {
70 if(gm2) free(gm2);
71 if(wt) free(wt);
72 if(mr) free(mr);
73 if(mg) free(mg);
74 if(mb) free(mb);
75 if(Qadd) free(Qadd);
76 throw FI_MSG_ERROR_MEMORY;
77 }
78 memset(gm2, 0, SIZE_3D * sizeof(float));
79 memset(wt, 0, SIZE_3D * sizeof(LONG));
80 memset(mr, 0, SIZE_3D * sizeof(LONG));
81 memset(mg, 0, SIZE_3D * sizeof(LONG));
82 memset(mb, 0, SIZE_3D * sizeof(LONG));
83 memset(Qadd, 0, sizeof(WORD) * width * height);
84 }
85
~WuQuantizer()86 WuQuantizer::~WuQuantizer() {
87 if(gm2) free(gm2);
88 if(wt) free(wt);
89 if(mr) free(mr);
90 if(mg) free(mg);
91 if(mb) free(mb);
92 if(Qadd) free(Qadd);
93 }
94
95
96 // Histogram is in elements 1..HISTSIZE along each axis,
97 // element 0 is for base or marginal value
98 // NB: these must start out 0!
99
100 // Build 3-D color histogram of counts, r/g/b, c^2
101 void
Hist3D(LONG * vwt,LONG * vmr,LONG * vmg,LONG * vmb,float * m2,int ReserveSize,RGBQUAD * ReservePalette)102 WuQuantizer::Hist3D(LONG *vwt, LONG *vmr, LONG *vmg, LONG *vmb, float *m2, int ReserveSize, RGBQUAD *ReservePalette) {
103 int ind = 0;
104 int inr, ing, inb, table[256];
105 int i;
106 unsigned y, x;
107
108 for(i = 0; i < 256; i++)
109 table[i] = i * i;
110
111 if (FreeImage_GetBPP(m_dib) == 24) {
112 for(y = 0; y < height; y++) {
113 BYTE *bits = FreeImage_GetScanLine(m_dib, y);
114
115 for(x = 0; x < width; x++) {
116 inr = (bits[FI_RGBA_RED] >> 3) + 1;
117 ing = (bits[FI_RGBA_GREEN] >> 3) + 1;
118 inb = (bits[FI_RGBA_BLUE] >> 3) + 1;
119 ind = INDEX(inr, ing, inb);
120 Qadd[y*width + x] = (WORD)ind;
121 // [inr][ing][inb]
122 vwt[ind]++;
123 vmr[ind] += bits[FI_RGBA_RED];
124 vmg[ind] += bits[FI_RGBA_GREEN];
125 vmb[ind] += bits[FI_RGBA_BLUE];
126 m2[ind] += (float)(table[bits[FI_RGBA_RED]] + table[bits[FI_RGBA_GREEN]] + table[bits[FI_RGBA_BLUE]]);
127 bits += 3;
128 }
129 }
130 } else {
131 for(y = 0; y < height; y++) {
132 BYTE *bits = FreeImage_GetScanLine(m_dib, y);
133
134 for(x = 0; x < width; x++) {
135 inr = (bits[FI_RGBA_RED] >> 3) + 1;
136 ing = (bits[FI_RGBA_GREEN] >> 3) + 1;
137 inb = (bits[FI_RGBA_BLUE] >> 3) + 1;
138 ind = INDEX(inr, ing, inb);
139 Qadd[y*width + x] = (WORD)ind;
140 // [inr][ing][inb]
141 vwt[ind]++;
142 vmr[ind] += bits[FI_RGBA_RED];
143 vmg[ind] += bits[FI_RGBA_GREEN];
144 vmb[ind] += bits[FI_RGBA_BLUE];
145 m2[ind] += (float)(table[bits[FI_RGBA_RED]] + table[bits[FI_RGBA_GREEN]] + table[bits[FI_RGBA_BLUE]]);
146 bits += 4;
147 }
148 }
149 }
150
151 if( ReserveSize > 0 ) {
152 int max = 0;
153 for(i = 0; i < SIZE_3D; i++) {
154 if( vwt[i] > max ) max = vwt[i];
155 }
156 max++;
157 for(i = 0; i < ReserveSize; i++) {
158 inr = (ReservePalette[i].rgbRed >> 3) + 1;
159 ing = (ReservePalette[i].rgbGreen >> 3) + 1;
160 inb = (ReservePalette[i].rgbBlue >> 3) + 1;
161 ind = INDEX(inr, ing, inb);
162 wt[ind] = max;
163 mr[ind] = max * ReservePalette[i].rgbRed;
164 mg[ind] = max * ReservePalette[i].rgbGreen;
165 mb[ind] = max * ReservePalette[i].rgbBlue;
166 gm2[ind] = (float)max * (float)(table[ReservePalette[i].rgbRed] + table[ReservePalette[i].rgbGreen] + table[ReservePalette[i].rgbBlue]);
167 }
168 }
169 }
170
171
172 // At conclusion of the histogram step, we can interpret
173 // wt[r][g][b] = sum over voxel of P(c)
174 // mr[r][g][b] = sum over voxel of r*P(c) , similarly for mg, mb
175 // m2[r][g][b] = sum over voxel of c^2*P(c)
176 // Actually each of these should be divided by 'ImageSize' to give the usual
177 // interpretation of P() as ranging from 0 to 1, but we needn't do that here.
178
179
180 // We now convert histogram into moments so that we can rapidly calculate
181 // the sums of the above quantities over any desired box.
182
183 // Compute cumulative moments
184 void
M3D(LONG * vwt,LONG * vmr,LONG * vmg,LONG * vmb,float * m2)185 WuQuantizer::M3D(LONG *vwt, LONG *vmr, LONG *vmg, LONG *vmb, float *m2) {
186 unsigned ind1, ind2;
187 BYTE i, r, g, b;
188 LONG line, line_r, line_g, line_b;
189 LONG area[33], area_r[33], area_g[33], area_b[33];
190 float line2, area2[33];
191
192 for(r = 1; r <= 32; r++) {
193 for(i = 0; i <= 32; i++) {
194 area2[i] = 0;
195 area[i] = area_r[i] = area_g[i] = area_b[i] = 0;
196 }
197 for(g = 1; g <= 32; g++) {
198 line2 = 0;
199 line = line_r = line_g = line_b = 0;
200 for(b = 1; b <= 32; b++) {
201 ind1 = INDEX(r, g, b); // [r][g][b]
202 line += vwt[ind1];
203 line_r += vmr[ind1];
204 line_g += vmg[ind1];
205 line_b += vmb[ind1];
206 line2 += m2[ind1];
207 area[b] += line;
208 area_r[b] += line_r;
209 area_g[b] += line_g;
210 area_b[b] += line_b;
211 area2[b] += line2;
212 ind2 = ind1 - 1089; // [r-1][g][b]
213 vwt[ind1] = vwt[ind2] + area[b];
214 vmr[ind1] = vmr[ind2] + area_r[b];
215 vmg[ind1] = vmg[ind2] + area_g[b];
216 vmb[ind1] = vmb[ind2] + area_b[b];
217 m2[ind1] = m2[ind2] + area2[b];
218 }
219 }
220 }
221 }
222
223 // Compute sum over a box of any given statistic
224 LONG
Vol(Box * cube,LONG * mmt)225 WuQuantizer::Vol( Box *cube, LONG *mmt ) {
226 return( mmt[INDEX(cube->r1, cube->g1, cube->b1)]
227 - mmt[INDEX(cube->r1, cube->g1, cube->b0)]
228 - mmt[INDEX(cube->r1, cube->g0, cube->b1)]
229 + mmt[INDEX(cube->r1, cube->g0, cube->b0)]
230 - mmt[INDEX(cube->r0, cube->g1, cube->b1)]
231 + mmt[INDEX(cube->r0, cube->g1, cube->b0)]
232 + mmt[INDEX(cube->r0, cube->g0, cube->b1)]
233 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
234 }
235
236 // The next two routines allow a slightly more efficient calculation
237 // of Vol() for a proposed subbox of a given box. The sum of Top()
238 // and Bottom() is the Vol() of a subbox split in the given direction
239 // and with the specified new upper bound.
240
241
242 // Compute part of Vol(cube, mmt) that doesn't depend on r1, g1, or b1
243 // (depending on dir)
244
245 LONG
Bottom(Box * cube,BYTE dir,LONG * mmt)246 WuQuantizer::Bottom(Box *cube, BYTE dir, LONG *mmt) {
247 switch(dir)
248 {
249 case FI_RGBA_RED:
250 return( - mmt[INDEX(cube->r0, cube->g1, cube->b1)]
251 + mmt[INDEX(cube->r0, cube->g1, cube->b0)]
252 + mmt[INDEX(cube->r0, cube->g0, cube->b1)]
253 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
254 break;
255 case FI_RGBA_GREEN:
256 return( - mmt[INDEX(cube->r1, cube->g0, cube->b1)]
257 + mmt[INDEX(cube->r1, cube->g0, cube->b0)]
258 + mmt[INDEX(cube->r0, cube->g0, cube->b1)]
259 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
260 break;
261 case FI_RGBA_BLUE:
262 return( - mmt[INDEX(cube->r1, cube->g1, cube->b0)]
263 + mmt[INDEX(cube->r1, cube->g0, cube->b0)]
264 + mmt[INDEX(cube->r0, cube->g1, cube->b0)]
265 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
266 break;
267 }
268
269 return 0;
270 }
271
272
273 // Compute remainder of Vol(cube, mmt), substituting pos for
274 // r1, g1, or b1 (depending on dir)
275
276 LONG
Top(Box * cube,BYTE dir,int pos,LONG * mmt)277 WuQuantizer::Top(Box *cube, BYTE dir, int pos, LONG *mmt) {
278 switch(dir)
279 {
280 case FI_RGBA_RED:
281 return( mmt[INDEX(pos, cube->g1, cube->b1)]
282 -mmt[INDEX(pos, cube->g1, cube->b0)]
283 -mmt[INDEX(pos, cube->g0, cube->b1)]
284 +mmt[INDEX(pos, cube->g0, cube->b0)] );
285 break;
286 case FI_RGBA_GREEN:
287 return( mmt[INDEX(cube->r1, pos, cube->b1)]
288 -mmt[INDEX(cube->r1, pos, cube->b0)]
289 -mmt[INDEX(cube->r0, pos, cube->b1)]
290 +mmt[INDEX(cube->r0, pos, cube->b0)] );
291 break;
292 case FI_RGBA_BLUE:
293 return( mmt[INDEX(cube->r1, cube->g1, pos)]
294 -mmt[INDEX(cube->r1, cube->g0, pos)]
295 -mmt[INDEX(cube->r0, cube->g1, pos)]
296 +mmt[INDEX(cube->r0, cube->g0, pos)] );
297 break;
298 }
299
300 return 0;
301 }
302
303 // Compute the weighted variance of a box
304 // NB: as with the raw statistics, this is really the variance * ImageSize
305
306 float
Var(Box * cube)307 WuQuantizer::Var(Box *cube) {
308 float dr = (float) Vol(cube, mr);
309 float dg = (float) Vol(cube, mg);
310 float db = (float) Vol(cube, mb);
311 float xx = gm2[INDEX(cube->r1, cube->g1, cube->b1)]
312 -gm2[INDEX(cube->r1, cube->g1, cube->b0)]
313 -gm2[INDEX(cube->r1, cube->g0, cube->b1)]
314 +gm2[INDEX(cube->r1, cube->g0, cube->b0)]
315 -gm2[INDEX(cube->r0, cube->g1, cube->b1)]
316 +gm2[INDEX(cube->r0, cube->g1, cube->b0)]
317 +gm2[INDEX(cube->r0, cube->g0, cube->b1)]
318 -gm2[INDEX(cube->r0, cube->g0, cube->b0)];
319
320 return (xx - (dr*dr+dg*dg+db*db)/(float)Vol(cube,wt));
321 }
322
323 // We want to minimize the sum of the variances of two subboxes.
324 // The sum(c^2) terms can be ignored since their sum over both subboxes
325 // is the same (the sum for the whole box) no matter where we split.
326 // The remaining terms have a minus sign in the variance formula,
327 // so we drop the minus sign and MAXIMIZE the sum of the two terms.
328
329 float
Maximize(Box * cube,BYTE dir,int first,int last,int * cut,LONG whole_r,LONG whole_g,LONG whole_b,LONG whole_w)330 WuQuantizer::Maximize(Box *cube, BYTE dir, int first, int last , int *cut, LONG whole_r, LONG whole_g, LONG whole_b, LONG whole_w) {
331 LONG half_r, half_g, half_b, half_w;
332 int i;
333 float temp;
334
335 LONG base_r = Bottom(cube, dir, mr);
336 LONG base_g = Bottom(cube, dir, mg);
337 LONG base_b = Bottom(cube, dir, mb);
338 LONG base_w = Bottom(cube, dir, wt);
339
340 float max = 0.0;
341
342 *cut = -1;
343
344 for (i = first; i < last; i++) {
345 half_r = base_r + Top(cube, dir, i, mr);
346 half_g = base_g + Top(cube, dir, i, mg);
347 half_b = base_b + Top(cube, dir, i, mb);
348 half_w = base_w + Top(cube, dir, i, wt);
349
350 // now half_x is sum over lower half of box, if split at i
351
352 if (half_w == 0) { // subbox could be empty of pixels!
353 continue; // never split into an empty box
354 } else {
355 temp = ((float)half_r*half_r + (float)half_g*half_g + (float)half_b*half_b)/half_w;
356 }
357
358 half_r = whole_r - half_r;
359 half_g = whole_g - half_g;
360 half_b = whole_b - half_b;
361 half_w = whole_w - half_w;
362
363 if (half_w == 0) { // subbox could be empty of pixels!
364 continue; // never split into an empty box
365 } else {
366 temp += ((float)half_r*half_r + (float)half_g*half_g + (float)half_b*half_b)/half_w;
367 }
368
369 if (temp > max) {
370 max=temp;
371 *cut=i;
372 }
373 }
374
375 return max;
376 }
377
378 bool
Cut(Box * set1,Box * set2)379 WuQuantizer::Cut(Box *set1, Box *set2) {
380 BYTE dir;
381 int cutr, cutg, cutb;
382
383 LONG whole_r = Vol(set1, mr);
384 LONG whole_g = Vol(set1, mg);
385 LONG whole_b = Vol(set1, mb);
386 LONG whole_w = Vol(set1, wt);
387
388 float maxr = Maximize(set1, FI_RGBA_RED, set1->r0+1, set1->r1, &cutr, whole_r, whole_g, whole_b, whole_w);
389 float maxg = Maximize(set1, FI_RGBA_GREEN, set1->g0+1, set1->g1, &cutg, whole_r, whole_g, whole_b, whole_w);
390 float maxb = Maximize(set1, FI_RGBA_BLUE, set1->b0+1, set1->b1, &cutb, whole_r, whole_g, whole_b, whole_w);
391
392 if ((maxr >= maxg) && (maxr >= maxb)) {
393 dir = FI_RGBA_RED;
394
395 if (cutr < 0) {
396 return false; // can't split the box
397 }
398 } else if ((maxg >= maxr) && (maxg>=maxb)) {
399 dir = FI_RGBA_GREEN;
400 } else {
401 dir = FI_RGBA_BLUE;
402 }
403
404 set2->r1 = set1->r1;
405 set2->g1 = set1->g1;
406 set2->b1 = set1->b1;
407
408 switch (dir) {
409 case FI_RGBA_RED:
410 set2->r0 = set1->r1 = cutr;
411 set2->g0 = set1->g0;
412 set2->b0 = set1->b0;
413 break;
414
415 case FI_RGBA_GREEN:
416 set2->g0 = set1->g1 = cutg;
417 set2->r0 = set1->r0;
418 set2->b0 = set1->b0;
419 break;
420
421 case FI_RGBA_BLUE:
422 set2->b0 = set1->b1 = cutb;
423 set2->r0 = set1->r0;
424 set2->g0 = set1->g0;
425 break;
426 }
427
428 set1->vol = (set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0);
429 set2->vol = (set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0);
430
431 return true;
432 }
433
434
435 void
Mark(Box * cube,int label,BYTE * tag)436 WuQuantizer::Mark(Box *cube, int label, BYTE *tag) {
437 for (int r = cube->r0 + 1; r <= cube->r1; r++) {
438 for (int g = cube->g0 + 1; g <= cube->g1; g++) {
439 for (int b = cube->b0 + 1; b <= cube->b1; b++) {
440 tag[INDEX(r, g, b)] = (BYTE)label;
441 }
442 }
443 }
444 }
445
446 // Wu Quantization algorithm
447 FIBITMAP *
Quantize(int PaletteSize,int ReserveSize,RGBQUAD * ReservePalette)448 WuQuantizer::Quantize(int PaletteSize, int ReserveSize, RGBQUAD *ReservePalette) {
449 BYTE *tag = NULL;
450
451 try {
452 Box cube[MAXCOLOR];
453 int next;
454 LONG i, weight;
455 int k;
456 float vv[MAXCOLOR], temp;
457
458 // Compute 3D histogram
459
460 Hist3D(wt, mr, mg, mb, gm2, ReserveSize, ReservePalette);
461
462 // Compute moments
463
464 M3D(wt, mr, mg, mb, gm2);
465
466 cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
467 cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
468 next = 0;
469
470 for (i = 1; i < PaletteSize; i++) {
471 if(Cut(&cube[next], &cube[i])) {
472 // volume test ensures we won't try to cut one-cell box
473 vv[next] = (cube[next].vol > 1) ? Var(&cube[next]) : 0;
474 vv[i] = (cube[i].vol > 1) ? Var(&cube[i]) : 0;
475 } else {
476 vv[next] = 0.0; // don't try to split this box again
477 i--; // didn't create box i
478 }
479
480 next = 0; temp = vv[0];
481
482 for (k = 1; k <= i; k++) {
483 if (vv[k] > temp) {
484 temp = vv[k]; next = k;
485 }
486 }
487
488 if (temp <= 0.0) {
489 PaletteSize = i + 1;
490
491 // Error: "Only got 'PaletteSize' boxes"
492
493 break;
494 }
495 }
496
497 // Partition done
498
499 // the space for array gm2 can be freed now
500
501 free(gm2);
502
503 gm2 = NULL;
504
505 // Allocate a new dib
506
507 FIBITMAP *new_dib = FreeImage_Allocate(width, height, 8);
508
509 if (new_dib == NULL) {
510 throw FI_MSG_ERROR_MEMORY;
511 }
512
513 // create an optimized palette
514
515 RGBQUAD *new_pal = FreeImage_GetPalette(new_dib);
516
517 tag = (BYTE*) malloc(SIZE_3D * sizeof(BYTE));
518 if (tag == NULL) {
519 throw FI_MSG_ERROR_MEMORY;
520 }
521 memset(tag, 0, SIZE_3D * sizeof(BYTE));
522
523 for (k = 0; k < PaletteSize ; k++) {
524 Mark(&cube[k], k, tag);
525 weight = Vol(&cube[k], wt);
526
527 if (weight) {
528 new_pal[k].rgbRed = (BYTE)(((float)Vol(&cube[k], mr) / (float)weight) + 0.5f);
529 new_pal[k].rgbGreen = (BYTE)(((float)Vol(&cube[k], mg) / (float)weight) + 0.5f);
530 new_pal[k].rgbBlue = (BYTE)(((float)Vol(&cube[k], mb) / (float)weight) + 0.5f);
531 } else {
532 // Error: bogus box 'k'
533
534 new_pal[k].rgbRed = new_pal[k].rgbGreen = new_pal[k].rgbBlue = 0;
535 }
536 }
537
538 int npitch = FreeImage_GetPitch(new_dib);
539
540 for (unsigned y = 0; y < height; y++) {
541 BYTE *new_bits = FreeImage_GetBits(new_dib) + (y * npitch);
542
543 for (unsigned x = 0; x < width; x++) {
544 new_bits[x] = tag[Qadd[y*width + x]];
545 }
546 }
547
548 // output 'new_pal' as color look-up table contents,
549 // 'new_bits' as the quantized image (array of table addresses).
550
551 free(tag);
552
553 return (FIBITMAP*) new_dib;
554 } catch(...) {
555 free(tag);
556 }
557
558 return NULL;
559 }
560