1// tQuantizeWu.cpp 
2// 
3// This module implements Wu quantization by Xiaolin Wu. The original header is included below. Modifications include: 
4// * Placing it in a namespace. 
5// * Consolidating the state parameters so that it is threadsafe (no global state). 
6// * Bridging to a standardized Tacent interface. 
7// * Convert to C++ syntax (original is pure C). 
8// * No exit or printf on error. 
9// 
10// The algrithm works well for larger numbers of colours (generally 128 to 256 or 255) but it can handle values as 
11// low as 2. 
12// 
13// Modifications Copyright (c) 2022-2024 Tristan Grimmer. 
14// Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby 
15// granted, provided that the above copyright notice and this permission notice appear in all copies. 
16// 
17// THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL 
18// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, 
19// INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
20// AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR 
21// PERFORMANCE OF THIS SOFTWARE. 
22// 
23// Here is the original code header from Xiaolin Wu: 
24// 
25// Having received many constructive comments and bug reports about my previous C implementation of my color quantizer 
26// (Graphics Gems vol. II, p. 126-133), I am posting the following second version of my program (hopefully 100% healthy) 
27// as a reply to all those who are interested in the problem. 
28// 
29// C Implementation of Wu's Color Quantizer (v. 2) (see Graphics Gems vol. II, pp. 126-133) 
30// Author: Xiaolin Wu, Dept. of Computer Science, Univ. of Western Ontario, London, Ontario N6A 5B7, wu@csd.uwo.ca 
31// Algorithm: Greedy orthogonal bipartition of RGB space for variance minimization aided by inclusion-exclusion tricks. 
32// For speed no nearest neighbor search is done. Slightly better performance can be expected by more sophisticated but 
33// more expensive versions. The author thanks Tom Lane at Tom_Lane@G.GP.CS.CMU.EDU for much of additional documentation 
34// and a cure to a previous bug. 
35// 
36// Free to distribute, comments and suggestions are appreciated. 
37#include <Math/tColour.h> 
38#include "Image/tQuantize.h" 
39namespace tImage { 
40 
41 
42namespace tQuantizeWu 
43{ 
44 // Constants. 
45 const int MaxColour = 256; // For 256 colours, fixed arrays need 8kb, plus space for the image. 
46 const int Red = 2; 
47 const int Green = 1; 
48 const int Blue = 0; 
49 
50 // We're putting all state on the stack (heap would work too). No globals. This allows the quantizer to be thread-safe. 
51 struct State 
52 { 
53 // Histogram is in elements 1..HISTSIZE along each axis, element 0 is for base or marginal value. These must start out 0. 
54 float m2[33][33][33]; 
55 int32 wt[33][33][33]; 
56 int32 mr[33][33][33]; 
57 int32 mg[33][33][33]; 
58 int32 mb[33][33][33]; 
59 uint8 *Ir, *Ig, *Ib; 
60 int size; // Image size (width*height). 
61 int K; // Colour look-up table size (palette size). 
62 uint16* Qadd; 
63 }; 
64 
65 struct Box 
66 { 
67 // X0 is min value, exclusive. X1 is max value, inclusive. 
68 int r0; int r1; int g0; int g1; int b0; int b1; 
69 int vol; 
70 }; 
71 
72 // Build 3-D color histogram of counts, r/g/b, c^2. At conclusion of the histogram step, we can interpret: 
73 // wt[r][g][b] = sum over voxel of P(c) 
74 // mr[r][g][b] = sum over voxel of r*P(c), similarly for mg, mb 
75 // m2[r][g][b] = sum over voxel of c^2*P(c) 
76 // Actually each of these should be divided by 'size' to give the usual interpretation of P() as ranging from 0 to 1 
77 // but we needn't do that here. 
78 void Hist3d(State&, int32* vwt, int32* vmr, int32* vmg, int32* vmb, float* m2); 
79 
80 // Compute cumulative moments. We now convert histogram into moments so that we can rapidly calculate the sums of 
81 // the above quantities over any desired box. 
82 void M3d(int32* vwt, int32* vmr, int32* vmg, int32* vmb, float* m2); 
83  
84 // Compute sum over a box of any given statistic. 
85 int32 Vol(Box* cube, int32 mmt[33][33][33]); 
86 
87 // The next two routines allow a slightly more efficient calculation of Vol() for a proposed subbox of a given box. 
88 // The sum of Top() and Bottom() is the Vol() of a subbox split in the given direction and with the specified new 
89 // upper bound. 
90 // 
91 // Compute part of Vol(cube, mmt) that doesn't depend on r1, g1, or b1 (depending on dir). 
92 int32 Bottom(Box* cube, uint8 dir, int32 mmt[33][33][33]); 
93 
94 // Compute remainder of Vol(cube, mmt), substituting pos for r1, g1, or b1 (depending on dir). 
95 int32 Top(Box* cube, uint8 dir, int pos, int32 mmt[33][33][33]); 
96 
97 // Compute the weighted variance of a box NB: as with the raw statistics, this is really the variance * size. 
98 float Var(State&, Box* cube); 
99 
100 // We want to minimize the sum of the variances of two subboxes. The sum(c^2) terms can be ignored since their sum 
101 // over both subboxes is the same (the sum for the whole box) no matter where we split. The remaining terms have a 
102 // minus sign in the variance formula, so we drop the minus sign and MAXIMIZE the sum of the two terms. 
103 float Maximize(State&, Box* cube, uint8 dir, int first, int last, int* cut, int32 whole_r, int32 whole_g, int32 whole_b, int32 whole_w); 
104 
105 int Cut(State&, Box* set1, Box* set2); 
106 void Mark(Box* cube, int label, uint8* tag); 
107 
108 // By the time this is called all parameters must be valid. It assumes the first 3 arguments are > 0 and the latter 3 non-null. 
109 void Quantize(int numColours, int width, int height, const tPixel3b* pixels, tColour3b* destPalette, uint8* destIndices); 
110} 
111 
112 
113void tQuantizeWu::Hist3d(State& state, int32* vwt, int32* vmr, int32* vmg, int32* vmb, float* m2) 
114{ 
115 int ind, r, g, b; 
116 int inr, ing, inb, table[256]; 
117 int32 i; 
118 
119 for (i = 0; i < 256; ++i) 
120 table[i] = i*i; 
121 
122 for (i = 0; i < state.size; ++i) 
123 { 
124 r = state.Ir[i]; g = state.Ig[i]; b = state.Ib[i]; 
125 inr = (r>>3)+1;  
126 ing = (g>>3)+1;  
127 inb = (b>>3)+1;  
128 state.Qadd[i] = ind = (inr<<10)+(inr<<6)+inr+(ing<<5)+ing+inb; 
129 // [inr][ing][inb] 
130 ++vwt[ind]; 
131 vmr[ind] += r; 
132 vmg[ind] += g; 
133 vmb[ind] += b; 
134 m2[ind] += (float)(table[r]+table[g]+table[b]); 
135 } 
136} 
137 
138 
139void tQuantizeWu::M3d(int32* vwt, int32* vmr, int32* vmg, int32* vmb, float* m2) 
140{ 
141 uint16 ind1, ind2; 
142 uint8 i, r, g, b; 
143 int32 line, line_r, line_g, line_b, area[33], area_r[33], area_g[33], area_b[33]; 
144 float line2, area2[33]; 
145 
146 for (r = 1; r <= 32; ++r) 
147 { 
148 for (i = 0; i <= 32; ++i) 
149 { 
150 area[i]=area_r[i]=area_g[i]=area_b[i]=0; 
151 area2[i] = 0.0f; 
152 } 
153 for (g = 1; g <= 32; ++g) 
154 { 
155 line = line_r = line_g = line_b = 0; 
156 line2 = 0.0f; 
157 for (b = 1; b <= 32; ++b) 
158 { 
159 ind1 = (r<<10) + (r<<6) + r + (g<<5) + g + b; // [r][g][b] 
160 line += vwt[ind1]; 
161 line_r += vmr[ind1];  
162 line_g += vmg[ind1];  
163 line_b += vmb[ind1]; 
164 line2 += m2[ind1]; 
165 area[b] += line; 
166 area_r[b] += line_r; 
167 area_g[b] += line_g; 
168 area_b[b] += line_b; 
169 area2[b] += line2; 
170 ind2 = ind1 - 1089; // [r-1][g][b] 
171 vwt[ind1] = vwt[ind2] + area[b]; 
172 vmr[ind1] = vmr[ind2] + area_r[b]; 
173 vmg[ind1] = vmg[ind2] + area_g[b]; 
174 vmb[ind1] = vmb[ind2] + area_b[b]; 
175 m2[ind1] = m2[ind2] + area2[b]; 
176 } 
177 } 
178 } 
179} 
180 
181 
182int32 tQuantizeWu::Vol(Box* cube, int32 mmt[33][33][33]) 
183{ 
184 return 
185 ( 
186 mmt[cube->r1][cube->g1][cube->b1] 
187 -mmt[cube->r1][cube->g1][cube->b0] 
188 -mmt[cube->r1][cube->g0][cube->b1] 
189 +mmt[cube->r1][cube->g0][cube->b0] 
190 -mmt[cube->r0][cube->g1][cube->b1] 
191 +mmt[cube->r0][cube->g1][cube->b0] 
192 +mmt[cube->r0][cube->g0][cube->b1] 
193 -mmt[cube->r0][cube->g0][cube->b0] 
194 ); 
195} 
196 
197 
198int32 tQuantizeWu::Bottom(Box* cube, uint8 dir, int32 mmt[33][33][33]) 
199{ 
200 switch(dir) 
201 { 
202 case Red: 
203 return 
204 ( 
205 -mmt[cube->r0][cube->g1][cube->b1] 
206 +mmt[cube->r0][cube->g1][cube->b0] 
207 +mmt[cube->r0][cube->g0][cube->b1] 
208 -mmt[cube->r0][cube->g0][cube->b0] 
209 ); 
210 
211 case Green: 
212 return 
213 ( 
214 -mmt[cube->r1][cube->g0][cube->b1] 
215 +mmt[cube->r1][cube->g0][cube->b0] 
216 +mmt[cube->r0][cube->g0][cube->b1] 
217 -mmt[cube->r0][cube->g0][cube->b0] 
218 ); 
219 
220 case Blue: 
221 return 
222 ( 
223 -mmt[cube->r1][cube->g1][cube->b0] 
224 +mmt[cube->r1][cube->g0][cube->b0] 
225 +mmt[cube->r0][cube->g1][cube->b0] 
226 -mmt[cube->r0][cube->g0][cube->b0] 
227 ); 
228 } 
229 
230 return 0; 
231} 
232 
233 
234int32 tQuantizeWu::Top(Box* cube, uint8 dir, int pos, int32 mmt[33][33][33]) 
235{ 
236 switch(dir) 
237 { 
238 case Red: 
239 return 
240 ( 
241 mmt[pos][cube->g1][cube->b1] 
242 -mmt[pos][cube->g1][cube->b0] 
243 -mmt[pos][cube->g0][cube->b1] 
244 +mmt[pos][cube->g0][cube->b0] 
245 ); 
246 
247 case Green: 
248 return 
249 ( 
250 mmt[cube->r1][pos][cube->b1] 
251 -mmt[cube->r1][pos][cube->b0] 
252 -mmt[cube->r0][pos][cube->b1] 
253 +mmt[cube->r0][pos][cube->b0] 
254 ); 
255 
256 case Blue: 
257 return 
258 ( 
259 mmt[cube->r1][cube->g1][pos] 
260 -mmt[cube->r1][cube->g0][pos] 
261 -mmt[cube->r0][cube->g1][pos] 
262 +mmt[cube->r0][cube->g0][pos] 
263 ); 
264 } 
265 return 0; 
266} 
267 
268 
269float tQuantizeWu::Var(State& state, Box* cube) 
270{ 
271 float dr, dg, db, xx; 
272 
273 dr = float( Vol(cube, mmt: state.mr) ); 
274 dg = float( Vol(cube, mmt: state.mg) ); 
275 db = float( Vol(cube, mmt: state.mb) ); 
276 xx = state.m2[cube->r1][cube->g1][cube->b1] 
277 -state.m2[cube->r1][cube->g1][cube->b0] 
278 -state.m2[cube->r1][cube->g0][cube->b1] 
279 +state.m2[cube->r1][cube->g0][cube->b0] 
280 -state.m2[cube->r0][cube->g1][cube->b1] 
281 +state.m2[cube->r0][cube->g1][cube->b0] 
282 +state.m2[cube->r0][cube->g0][cube->b1] 
283 -state.m2[cube->r0][cube->g0][cube->b0]; 
284 
285 return (xx - (dr*dr+dg*dg+db*db)/(float)Vol(cube,mmt: state.wt)); 
286} 
287 
288 
289float tQuantizeWu::Maximize(State& state, Box* cube, uint8 dir, int first, int last, int* cut, int32 whole_r, int32 whole_g, int32 whole_b, int32 whole_w) 
290{ 
291 int32 half_r, half_g, half_b, half_w; 
292 int32 base_r, base_g, base_b, base_w; 
293 int i; 
294 float temp, max; 
295 
296 base_r = Bottom(cube, dir, mmt: state.mr); 
297 base_g = Bottom(cube, dir, mmt: state.mg); 
298 base_b = Bottom(cube, dir, mmt: state.mb); 
299 base_w = Bottom(cube, dir, mmt: state.wt); 
300 max = 0.0; 
301 *cut = -1; 
302 for (i = first; i < last; ++i) 
303 { 
304 half_r = base_r + Top(cube, dir, pos: i, mmt: state.mr); 
305 half_g = base_g + Top(cube, dir, pos: i, mmt: state.mg); 
306 half_b = base_b + Top(cube, dir, pos: i, mmt: state.mb); 
307 half_w = base_w + Top(cube, dir, pos: i, mmt: state.wt); 
308 
309 // Now half_x is sum over lower half of box, if split at i. 
310 // Subbox could be empty of pixels. Never split into an empty box. 
311 if (half_w == 0) 
312 continue; 
313 else 
314 temp = ((float)half_r*half_r + (float)half_g*half_g + (float)half_b*half_b)/half_w; 
315 
316 half_r = whole_r - half_r; 
317 half_g = whole_g - half_g; 
318 half_b = whole_b - half_b; 
319 half_w = whole_w - half_w; 
320 
321 // Subbox could be empty of pixels. Never split into an empty box. 
322 if (half_w == 0) 
323 continue; 
324 else 
325 temp += ((float)half_r*half_r + (float)half_g*half_g + (float)half_b*half_b)/half_w; 
326 
327 if (temp > max) 
328 { 
329 max=temp; 
330 *cut=i; 
331 } 
332 } 
333 
334 return max; 
335} 
336 
337 
338int tQuantizeWu::Cut(State& state, Box* set1, Box* set2) 
339{ 
340 uint8 dir; 
341 int32 cutr, cutg, cutb; 
342 float maxr, maxg, maxb; 
343 int32 whole_r, whole_g, whole_b, whole_w; 
344 
345 whole_r = Vol(cube: set1, mmt: state.mr); 
346 whole_g = Vol(cube: set1, mmt: state.mg); 
347 whole_b = Vol(cube: set1, mmt: state.mb); 
348 whole_w = Vol(cube: set1, mmt: state.wt); 
349 
350 maxr = Maximize(state, cube: set1, dir: Red, first: set1->r0+1, last: set1->r1, cut: &cutr, whole_r, whole_g, whole_b, whole_w); 
351 maxg = Maximize(state, cube: set1, dir: Green, first: set1->g0+1, last: set1->g1, cut: &cutg, whole_r, whole_g, whole_b, whole_w); 
352 maxb = Maximize(state, cube: set1, dir: Blue, first: set1->b0+1, last: set1->b1, cut: &cutb, whole_r, whole_g, whole_b, whole_w); 
353 
354 if ((maxr >= maxg) && (maxr >= maxb)) 
355 { 
356 dir = Red; 
357 
358 // Can't split the box. 
359 if (cutr < 0) 
360 return 0; 
361 } 
362 else if ((maxg >= maxr) && (maxg >= maxb)) 
363 dir = Green; 
364 else 
365 dir = Blue; 
366 
367 set2->r1 = set1->r1; 
368 set2->g1 = set1->g1; 
369 set2->b1 = set1->b1; 
370 
371 switch (dir) 
372 { 
373 case Red: 
374 set2->r0 = set1->r1 = cutr; 
375 set2->g0 = set1->g0; 
376 set2->b0 = set1->b0; 
377 break; 
378 
379 case Green: 
380 set2->g0 = set1->g1 = cutg; 
381 set2->r0 = set1->r0; 
382 set2->b0 = set1->b0; 
383 break; 
384 
385 case Blue: 
386 set2->b0 = set1->b1 = cutb; 
387 set2->r0 = set1->r0; 
388 set2->g0 = set1->g0; 
389 break; 
390 } 
391 
392 set1->vol=(set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0); 
393 set2->vol=(set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0); 
394 return 1; 
395} 
396 
397 
398void tQuantizeWu::Mark(Box* cube, int label, uint8* tag) 
399{ 
400 int r, g, b; 
401 for(r=cube->r0+1; r<=cube->r1; ++r) 
402 for(g=cube->g0+1; g<=cube->g1; ++g) 
403 for(b=cube->b0+1; b<=cube->b1; ++b) 
404 tag[(r<<10) + (r<<6) + r + (g<<5) + g + b] = label; 
405} 
406 
407 
408void tQuantizeWu::Quantize(int numColours, int width, int height, const tPixel3b* pixels, tColour3b* destPalette, uint8* destIndices) 
409{ 
410 Box cube[MaxColour]; 
411 uint8 lut_r[MaxColour], lut_g[MaxColour], lut_b[MaxColour]; 
412 int next; 
413 int32 i, weight; 
414 int k; 
415 float vv[MaxColour], temp; 
416 int numPixels = width*height; 
417 
418 // The arrays in the state are required to be initialized to 0 (and 0.0f for the float array). 
419 State state; 
420 tStd::tMemset(dest: &state, val: 0, numBytes: sizeof(State)); 
421 
422 state.K = numColours; 
423 state.size = numPixels; 
424 
425 // This is where the indices will be stored. 
426 state.Qadd = new uint16[numPixels]; 
427 
428 // Input R,G,B components into Ir, Ig, Ib. 
429 state.Ir = new uint8[numPixels]; 
430 state.Ig = new uint8[numPixels]; 
431 state.Ib = new uint8[numPixels]; 
432 for (int y = 0; y < height; y++) 
433 { 
434 for (int x = 0; x < width; x++) 
435 { 
436 int index = x + y*width; 
437 const tPixel3b& pixel = pixels[index]; 
438 state.Ir[index] = pixel.R; 
439 state.Ig[index] = pixel.G; 
440 state.Ib[index] = pixel.B; 
441 } 
442 } 
443 
444 Hist3d(state, vwt: (int32*)state.wt, vmr: (int32*)state.mr, vmg: (int32*)state.mg, vmb: (int32*)state.mb, m2: (float*)state.m2); 
445  
446 // These won't be accessed again. 
447 delete[] state.Ib; 
448 delete[] state.Ig; 
449 delete[] state.Ir; 
450 
451 M3d(vwt: (int32*)state.wt, vmr: (int32*)state.mr, vmg: (int32*)state.mg, vmb: (int32*)state.mb, m2: (float*)state.m2); 
452 cube[0].r0 = cube[0].g0 = cube[0].b0 = 0; 
453 cube[0].r1 = cube[0].g1 = cube[0].b1 = 32; 
454 next = 0; 
455 
456 for (i = 1; i < state.K; ++i) 
457 { 
458 if (Cut(state, set1: &cube[next], set2: &cube[i])) 
459 { 
460 // Volume test ensures we won't try to cut one-cell box. 
461 vv[next] = (cube[next].vol>1) ? Var(state, cube: &cube[next]) : 0.0; 
462 vv[i] = (cube[i].vol>1) ? Var(state, cube: &cube[i]) : 0.0; 
463 } 
464 else 
465 { 
466 vv[next] = 0.0; // Don't try to split this box again. 
467 i--; // Didn't create box i. 
468 } 
469 
470 next = 0; temp = vv[0]; 
471 for (k = 1; k <= i; ++k) 
472 { 
473 if (vv[k] > temp) 
474 { 
475 temp = vv[k]; 
476 next = k; 
477 } 
478 } 
479 
480 if (temp <= 0.0) 
481 { 
482 state.K = i+1; 
483 break; 
484 } 
485 } 
486 
487 // The m2 float array will not be accessed further now. 
488 uint8 tag[33*33*33]; 
489 for (k = 0; k < state.K; ++k) 
490 { 
491 Mark(cube: &cube[k], label: k, tag); 
492 weight = Vol(cube: &cube[k], mmt: state.wt); 
493 if (weight) 
494 { 
495 lut_r[k] = Vol(cube: &cube[k], mmt: state.mr) / weight; 
496 lut_g[k] = Vol(cube: &cube[k], mmt: state.mg) / weight; 
497 lut_b[k] = Vol(cube: &cube[k], mmt: state.mb) / weight; 
498 } 
499 else 
500 { 
501 // k is a bogus box. 
502 lut_r[k] = lut_g[k] = lut_b[k] = 0;  
503 } 
504 } 
505 
506 for (i = 0; i < state.size; ++i) 
507 state.Qadd[i] = tag[state.Qadd[i]]; 
508 
509 // Populate the palette. lut_r, lut_g, and lut_b contain the lookup table colours. 
510 tAssert(state.K <= numColours); 
511 for (int ind = 0; ind < state.K; ind++) 
512 { 
513 destPalette[ind].Set(r: lut_r[ind], g: lut_g[ind], b: lut_b[ind]); 
514 } 
515 
516 // Copy the indices into the supplied dest array. These are stored in Qadd. 
517 for (int ind = 0; ind < numPixels; ind++) 
518 { 
519 tAssert( state.Qadd[ind] <= 255 ); 
520 destIndices[ind] = uint8(state.Qadd[ind]); 
521 } 
522  
523 delete[] state.Qadd; 
524} 
525 
526 
527// 
528// The functions below make up the external interface. 
529// 
530 
531 
532bool tQuantizeWu::QuantizeImage 
533( 
534 int numColours, int width, int height, const tPixel3b* pixels, tColour3b* destPalette, uint8* destIndices, 
535 bool checkExact 
536) 
537{ 
538 if ((numColours < 2) || (numColours > 256) || (width <= 0) || (height <= 0) || !pixels || !destPalette || !destIndices) 
539 return false; 
540 
541 if (checkExact) 
542 { 
543 bool success = tQuantize::QuantizeImageExact(numColours, width, height, pixels, destPalette, destIndices); 
544 if (success) 
545 return true; 
546 } 
547 
548 Quantize(numColours, width, height, pixels, destPalette, destIndices); 
549 return true; 
550} 
551 
552 
553bool tQuantizeWu::QuantizeImage 
554( 
555 int numColours, int width, int height, const tPixel4b* pixels, tColour3b* destPalette, uint8* destIndices, 
556 bool checkExact 
557) 
558{ 
559 if ((numColours < 2) || (numColours > 256) || (width <= 0) || (height <= 0) || !pixels || !destPalette || !destIndices) 
560 return false; 
561 
562 tPixel3b* pixels3 = new tPixel3b[width*height]; 
563 for (int p = 0; p < width*height; p++) 
564 pixels3[p].Set( r: pixels[p].R, g: pixels[p].G, b: pixels[p].B ); 
565 
566 bool success = QuantizeImage(numColours, width, height, pixels: pixels3, destPalette, destIndices, checkExact); 
567 delete[] pixels3; 
568 return success; 
569} 
570 
571 
572} 
573