1 /********************************************************************************
2 * *
3 * W U C o l o r Q u a n t i z a t i o n *
4 * *
5 *********************************************************************************
6 * Copyright (C) 2004,2021 by Jeroen van der Zijp. All Rights Reserved. *
7 *********************************************************************************
8 * This library is free software; you can redistribute it and/or modify *
9 * it under the terms of the GNU Lesser General Public License as published by *
10 * the Free Software Foundation; either version 3 of the License, or *
11 * (at your option) any later version. *
12 * *
13 * This library is distributed in the hope that it will be useful, *
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16 * GNU Lesser General Public License for more details. *
17 * *
18 * You should have received a copy of the GNU Lesser General Public License *
19 * along with this program. If not, see <http://www.gnu.org/licenses/> *
20 ********************************************************************************/
21 #include "xincs.h"
22 #include "fxver.h"
23 #include "fxdefs.h"
24 #include "fxmath.h"
25
26 /*
27 Notes:
28
29 - This code is due to: Xiaolin Wu, Dept. of Computer Science, Univ. of
30 Western Ontario, London, Ontario N6A 5B7 (wu@csd.uwo.ca).
31 Original algorithm can be found in Graphics Gems vol. II, pp. 126-133.
32
33 - Algorithm: Greedy orthogonal bipartition of RGB space for variance
34 minimization aided by inclusion-exclusion tricks.
35 For speed no nearest neighbor search is done. Slightly
36 better performance can be expected by more sophisticated
37 but more expensive versions.
38
39 - Modified by Jeroen for FOX; don't blame the original author if I
40 broke it.
41
42 */
43
44 #define MAXCOLOR 256
45 #define RED 2
46 #define GREEN 1
47 #define BLUE 0
48
49
50 using namespace FX;
51
52 /*******************************************************************************/
53
54 namespace FX {
55
56
57 extern FXbool fxwuquantize(FXuchar* dst,const FXColor* src,FXColor* colormap,FXint& actualcolors,FXint w,FXint h,FXint maxcolors);
58
59
60 // Sub bpx
61 struct box {
62 FXint r0; // Min value, exclusive
63 FXint r1; // Max value, inclusive
64 FXint g0;
65 FXint g1;
66 FXint b0;
67 FXint b1;
68 FXint vol; // Volume
69 };
70
71
72 // To pass around
73 struct WU {
74 FXfloat m2[33][33][33]; // Histogram is in elements 1..HISTSIZE along each axis,
75 FXint wt[33][33][33]; // element 0 is for base or marginal value
76 FXint mr[33][33][33]; // NB: these must start out 0!
77 FXint mg[33][33][33];
78 FXint mb[33][33][33];
79 };
80
81
82
83 // At conclusion of the histogram step, we can interpret
84 // wt[r][g][b] = sum over voxel of P(c)
85 // mr[r][g][b] = sum over voxel of r*P(c) , similarly for mg, mb
86 // m2[r][g][b] = sum over voxel of c^2*P(c)
87 // Actually each of these should be divided by 'size' to give the usual
88 // interpretation of P() as ranging from 0 to 1, but we needn't do that here.
89
90 // Build 3-D color histogram of counts, r/g/b, c^2
histogram(WU & wu,const FXColor * data,FXint size)91 static void histogram(WU& wu,const FXColor *data,FXint size){
92 FXint r,g,b,inr,ing,inb,i;
93
94 // Clear counters
95 memset(&wu,0,sizeof(wu));
96
97 // Build histogram
98 for(i=0; i<size; ++i){
99 r=((const FXuchar*)(data+i))[2];
100 g=((const FXuchar*)(data+i))[1];
101 b=((const FXuchar*)(data+i))[0];
102 inr=(r>>3)+1;
103 ing=(g>>3)+1;
104 inb=(b>>3)+1;
105 wu.wt[inr][ing][inb]+=1;
106 wu.mr[inr][ing][inb]+=r;
107 wu.mg[inr][ing][inb]+=g;
108 wu.mb[inr][ing][inb]+=b;
109 wu.m2[inr][ing][inb]+=(FXfloat)(r*r+g*g+b*b);
110 }
111 }
112
113
114 // Compute cumulative moments
moments(WU & wu)115 static void moments(WU& wu){
116 FXint linet,liner,lineg,lineb,i,r,g,b;
117 FXint areat[33],arear[33],areag[33],areab[33];
118 FXfloat line2,area2[33];
119 for(r=1; r<=32; ++r){
120 for(i=0; i<=32; ++i){
121 areat[i]=0;
122 arear[i]=0;
123 areag[i]=0;
124 areab[i]=0;
125 area2[i]=0.0f;
126 }
127 for(g=1; g<=32; ++g){
128 linet=0;
129 liner=0;
130 lineg=0;
131 lineb=0;
132 line2=0.0f;
133 for(b=1; b<=32; ++b){
134 linet+=wu.wt[r][g][b];
135 liner+=wu.mr[r][g][b];
136 lineg+=wu.mg[r][g][b];
137 lineb+=wu.mb[r][g][b];
138 line2+=wu.m2[r][g][b];
139 areat[b]+=linet;
140 arear[b]+=liner;
141 areag[b]+=lineg;
142 areab[b]+=lineb;
143 area2[b]+=line2;
144 wu.wt[r][g][b]=wu.wt[r-1][g][b]+areat[b];
145 wu.mr[r][g][b]=wu.mr[r-1][g][b]+arear[b];
146 wu.mg[r][g][b]=wu.mg[r-1][g][b]+areag[b];
147 wu.mb[r][g][b]=wu.mb[r-1][g][b]+areab[b];
148 wu.m2[r][g][b]=wu.m2[r-1][g][b]+area2[b];
149 }
150 }
151 }
152 }
153
154
155 // Compute sum over a box of any given statistic
volume(box & cube,FXint mmt[33][33][33])156 static int volume(box& cube,FXint mmt[33][33][33]){
157 return mmt[cube.r1][cube.g1][cube.b1]
158 -mmt[cube.r1][cube.g1][cube.b0]
159 -mmt[cube.r1][cube.g0][cube.b1]
160 +mmt[cube.r1][cube.g0][cube.b0]
161 -mmt[cube.r0][cube.g1][cube.b1]
162 +mmt[cube.r0][cube.g1][cube.b0]
163 +mmt[cube.r0][cube.g0][cube.b1]
164 -mmt[cube.r0][cube.g0][cube.b0];
165 }
166
167
168 // The next two routines allow a slightly more efficient calculation
169 // of Vol() for a proposed subbox of a given box. The sum of Top()
170 // and Bottom() is the Vol() of a subbox split in the given direction
171 // and with the specified new upper bound.
172
173 // Compute part of Vol(cube, mmt) that doesn't depend
174 // on r1, g1, or b1 (depending on dir)
bottom(box & cube,FXuchar dir,FXint mmt[33][33][33])175 static FXint bottom(box& cube,FXuchar dir,FXint mmt[33][33][33]){
176 FXint result=0;
177 switch(dir){
178 case RED:
179 result= -mmt[cube.r0][cube.g1][cube.b1]
180 +mmt[cube.r0][cube.g1][cube.b0]
181 +mmt[cube.r0][cube.g0][cube.b1]
182 -mmt[cube.r0][cube.g0][cube.b0];
183 break;
184 case GREEN:
185 result= -mmt[cube.r1][cube.g0][cube.b1]
186 +mmt[cube.r1][cube.g0][cube.b0]
187 +mmt[cube.r0][cube.g0][cube.b1]
188 -mmt[cube.r0][cube.g0][cube.b0];
189 break;
190 case BLUE:
191 result= -mmt[cube.r1][cube.g1][cube.b0]
192 +mmt[cube.r1][cube.g0][cube.b0]
193 +mmt[cube.r0][cube.g1][cube.b0]
194 -mmt[cube.r0][cube.g0][cube.b0];
195 break;
196 }
197 return result;
198 }
199
200
201 // Compute remainder of Vol(cube, mmt), substituting pos
202 // for r1, g1, or b1 (depending on dir)
top(box & cube,FXuchar dir,FXint pos,FXint mmt[33][33][33])203 static FXint top(box& cube,FXuchar dir,FXint pos,FXint mmt[33][33][33]){
204 FXint result=0;
205 switch(dir){
206 case RED:
207 result= mmt[pos][cube.g1][cube.b1]
208 -mmt[pos][cube.g1][cube.b0]
209 -mmt[pos][cube.g0][cube.b1]
210 +mmt[pos][cube.g0][cube.b0];
211 break;
212 case GREEN:
213 result= mmt[cube.r1][pos][cube.b1]
214 -mmt[cube.r1][pos][cube.b0]
215 -mmt[cube.r0][pos][cube.b1]
216 +mmt[cube.r0][pos][cube.b0];
217 break;
218 case BLUE:
219 result= mmt[cube.r1][cube.g1][pos]
220 -mmt[cube.r1][cube.g0][pos]
221 -mmt[cube.r0][cube.g1][pos]
222 +mmt[cube.r0][cube.g0][pos];
223 break;
224 }
225 return result;
226 }
227
228
229 // Compute the weighted variance of a box
230 // NB: as with the raw statistics, this is really the variance * size
variance(WU & wu,box & cube)231 static FXfloat variance(WU& wu,box& cube){
232 FXfloat dr,dg,db,xx;
233
234 dr = (FXfloat)volume(cube,wu.mr);
235 dg = (FXfloat)volume(cube,wu.mg);
236 db = (FXfloat)volume(cube,wu.mb);
237
238 xx = wu.m2[cube.r1][cube.g1][cube.b1]
239 -wu.m2[cube.r1][cube.g1][cube.b0]
240 -wu.m2[cube.r1][cube.g0][cube.b1]
241 +wu.m2[cube.r1][cube.g0][cube.b0]
242 -wu.m2[cube.r0][cube.g1][cube.b1]
243 +wu.m2[cube.r0][cube.g1][cube.b0]
244 +wu.m2[cube.r0][cube.g0][cube.b1]
245 -wu.m2[cube.r0][cube.g0][cube.b0];
246
247 return xx-(dr*dr+dg*dg+db*db)/(FXfloat)volume(cube,wu.wt);
248 }
249
250
251 // We want to minimize the sum of the variances of two subboxes.
252 // The sum(c^2) terms can be ignored since their sum over both subboxes
253 // is the same (the sum for the whole box) no matter where we split.
254 // The remaining terms have a minus sign in the variance formula,
255 // so we drop the minus sign and MAXIMIZE the sum of the two terms.
maximize(WU & wu,box & cube,FXuchar dir,FXint first,FXint last,FXint * cut,FXint whole_r,FXint whole_g,FXint whole_b,FXint whole_w)256 static FXfloat maximize(WU& wu,box& cube,FXuchar dir,FXint first,FXint last,FXint *cut,FXint whole_r,FXint whole_g,FXint whole_b,FXint whole_w){
257 FXint half_r,half_g,half_b,half_w,base_r,base_g,base_b,base_w,i;
258 FXfloat temp,max;
259
260 base_r=bottom(cube,dir,wu.mr);
261 base_g=bottom(cube,dir,wu.mg);
262 base_b=bottom(cube,dir,wu.mb);
263 base_w=bottom(cube,dir,wu.wt);
264
265 max=0.0f;
266 *cut=-1;
267 for(i=first; i<last; ++i){
268
269 half_r=base_r+top(cube,dir,i,wu.mr);
270 half_g=base_g+top(cube,dir,i,wu.mg);
271 half_b=base_b+top(cube,dir,i,wu.mb);
272 half_w=base_w+top(cube,dir,i,wu.wt);
273
274 // Now half_x is sum over lower half of box, if split at i
275
276 // Subbox could be empty of pixels; never split into an empty box
277 if(half_w==0) continue;
278
279 temp=((FXfloat)half_r*half_r + (FXfloat)half_g*half_g + (FXfloat)half_b*half_b)/half_w;
280
281 half_r=whole_r-half_r;
282 half_g=whole_g-half_g;
283 half_b=whole_b-half_b;
284 half_w=whole_w-half_w;
285
286 // Subbox could be empty of pixels; never split into an empty box
287 if(half_w==0) continue;
288
289 temp += ((FXfloat)half_r*half_r + (FXfloat)half_g*half_g + (FXfloat)half_b*half_b)/half_w;
290
291 if(temp>max){ max=temp; *cut=i; }
292 }
293 return max;
294 }
295
296
297 // Find best split
cut(WU & wu,box & set1,box & set2)298 static FXint cut(WU& wu,box& set1,box& set2){
299 FXint cutr, cutg, cutb;
300 FXfloat maxr, maxg, maxb;
301 FXint whole_r, whole_g, whole_b, whole_w;
302 FXuchar dir;
303
304 // Totals
305 whole_r=volume(set1,wu.mr);
306 whole_g=volume(set1,wu.mg);
307 whole_b=volume(set1,wu.mb);
308 whole_w=volume(set1,wu.wt);
309
310 // Find most beneficial split direction
311 maxr=maximize(wu,set1, RED,set1.r0+1,set1.r1,&cutr,whole_r,whole_g,whole_b,whole_w);
312 maxg=maximize(wu,set1,GREEN,set1.g0+1,set1.g1,&cutg,whole_r,whole_g,whole_b,whole_w);
313 maxb=maximize(wu,set1, BLUE,set1.b0+1,set1.b1,&cutb,whole_r,whole_g,whole_b,whole_w);
314
315 // Direction of split?
316 if((maxr>=maxg) && (maxr>=maxb)){
317 if(cutr<0) return 0; // Can't split the box
318 dir=RED;
319 }
320 else if((maxg>=maxr) && (maxg>=maxb)){
321 dir=GREEN;
322 }
323 else{
324 dir=BLUE;
325 }
326
327 set2.r1=set1.r1;
328 set2.g1=set1.g1;
329 set2.b1=set1.b1;
330
331 switch(dir){
332 case RED:
333 set2.r0=set1.r1=cutr;
334 set2.g0=set1.g0;
335 set2.b0=set1.b0;
336 break;
337 case GREEN:
338 set2.g0=set1.g1=cutg;
339 set2.r0=set1.r0;
340 set2.b0=set1.b0;
341 break;
342 case BLUE:
343 set2.b0=set1.b1=cutb;
344 set2.r0=set1.r0;
345 set2.g0=set1.g0;
346 break;
347 }
348
349 set1.vol=(set1.r1-set1.r0)*(set1.g1-set1.g0)*(set1.b1-set1.b0);
350 set2.vol=(set2.r1-set2.r0)*(set2.g1-set2.g0)*(set2.b1-set2.b0);
351 return 1;
352 }
353
354
355 // Each entry in box maps to label
mark(box & cube,FXint label,FXuchar map[33][33][33])356 static void mark(box& cube,FXint label,FXuchar map[33][33][33]){
357 FXint r,g,b;
358 for(r=cube.r0+1; r<=cube.r1; ++r){
359 for(g=cube.g0+1; g<=cube.g1; ++g){
360 for(b=cube.b0+1; b<=cube.b1; ++b){
361 map[r][g][b]=label;
362 }
363 }
364 }
365 }
366
367
368 // Wu's quantization method based on recursive partitioning
fxwuquantize(FXuchar * dst,const FXColor * src,FXColor * colormap,FXint & actualcolors,FXint w,FXint h,FXint maxcolors)369 FXbool fxwuquantize(FXuchar* dst,const FXColor* src,FXColor* colormap,FXint& actualcolors,FXint w,FXint h,FXint maxcolors){
370 FXint i,k,weight,next,size,r,g,b;
371 FXfloat temp;
372 FXuchar map[33][33][33];
373 FXfloat vv[MAXCOLOR];
374 box cube[MAXCOLOR];
375 WU wu;
376
377 // Size of image
378 size=w*h;
379
380 // Compute histogram
381 histogram(wu,src,size);
382
383 // Compute moments
384 moments(wu);
385
386 // Recursively split boxes
387 next=0;
388 cube[0].r0=cube[0].g0=cube[0].b0=0;
389 cube[0].r1=cube[0].g1=cube[0].b1=32;
390 for(i=1; i<maxcolors; ++i){
391 if(cut(wu,cube[next],cube[i])){
392 vv[next]=(cube[next].vol>1)?variance(wu,cube[next]):0.0f; // Volume test ensures we won't try to cut one-cell box
393 vv[i]=(cube[i].vol>1)?variance(wu,cube[i]):0.0f;
394 }
395 else{
396 vv[next]=0.0f; // Don't try to split this box again
397 i--; // Didn't create box i
398 }
399 next=0;
400 temp=vv[0];
401 for(k=1; k<=i; ++k){
402 if(vv[k]>temp){
403 temp=vv[k];
404 next=k;
405 }
406 }
407 if(temp<=0.0f){
408 maxcolors=i+1;
409 break;
410 }
411 }
412
413 // Construct colormap
414 for(k=0; k<maxcolors; ++k){
415 mark(cube[k],k,map);
416 weight=volume(cube[k],wu.wt);
417 if(weight){
418 ((FXuchar*)(colormap+k))[3]=255;
419 ((FXuchar*)(colormap+k))[2]=volume(cube[k],wu.mr)/weight;
420 ((FXuchar*)(colormap+k))[1]=volume(cube[k],wu.mg)/weight;
421 ((FXuchar*)(colormap+k))[0]=volume(cube[k],wu.mb)/weight;
422 }
423 else{
424 ((FXuchar*)(colormap+k))[0]=0;
425 ((FXuchar*)(colormap+k))[1]=0;
426 ((FXuchar*)(colormap+k))[2]=0;
427 ((FXuchar*)(colormap+k))[3]=0;
428 }
429 }
430
431 // Quantize image
432 for(i=0; i<size; ++i){
433 r=((const FXuchar*)(src+i))[2];
434 g=((const FXuchar*)(src+i))[1];
435 b=((const FXuchar*)(src+i))[0];
436 dst[i]=map[(r>>3)+1][(g>>3)+1][(b>>3)+1];
437 }
438
439 // Return actual number of colors
440 actualcolors=maxcolors;
441
442 return true;
443 }
444
445 }
446