|
/Users/lyon/j4p/src/ip/color/Wu.java
|
1 /*
2 Master's Project step1
3 Color quantization using Wu's method
4 "Efficient Statistical Computations for
5 Optimal Color Quantization"
6 in Graphics Gem Vol II, edited by James Arvo, AP.
7 pp 116-125
8 Based on C code available from:
9 <http://www.cis.ohio-state.edu/~parent/classes/781/ip.graphics/GG2>
10
11 Hak Kywn Roh
12
13 */
14
15 package ip.color;
16
17
18 class Box {
19 public int r0; /* min value, exclusive */
20 public int r1; /* max value, inclusive */
21 public int g0;
22 public int g1;
23 public int b0;
24 public int b1;
25 public int vol;
26
27 public Box() {
28 }
29
30
31 }
32
33 public class Wu {
34 static final int MAXCOLOR = 256;
35 static final int RED = 2;
36 static final int GREEN = 1;
37 static final int BLUE = 0;
38
39 int size, K;
40 int qs = 33; //quant size
41 float m2[][][] = new float[qs][qs][qs];
42 int wt[][][] = new int[qs][qs][qs];
43 int mr[][][] = new int[qs][qs][qs];
44 int mg[][][] = new int[qs][qs][qs];
45 int mb[][][] = new int[qs][qs][qs];
46 int qadd[][];
47 short ir[][], ig[][], ib[][];
48 int height, width;
49
50
51 public void wuQuantization(short r[][],
52 short g[][],
53 short b[][],
54 int w,
55 int h,
56 int nK) {
57 width = w;
58 height = h;
59 ir = r;
60 ig = g;
61 ib = b;
62 size = width * height;
63 K = nK;
64
65 int next, i, k, x, y, z;
66 Box[] cube = new Box[MAXCOLOR];
67 int tag[][][] = new int[qs][qs][qs];
68 short lutr[] = new short[MAXCOLOR];
69 short lutg[] = new short[MAXCOLOR];
70 short lutb[] = new short[MAXCOLOR];
71 float vv[] = new float[MAXCOLOR];
72 float temp;
73 long weight;
74
75 Hist3d(wt, mr, mg, mb, m2);
76 System.out.println("Histogram done...");
77
78 M3d(wt, mr, mg, mb, m2);
79 System.out.println("Moments done...");
80
81 for (z = 0; z < MAXCOLOR; z++)
82 cube[z] = new Box();
83
84 cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
85 cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
86
87 next = 0;
88
89 for (i = 1; i < K; i++) {
90 if (Cut(cube[next], cube[i]) != -1) {
91 if (cube[next].vol > 1)
92 vv[next] = Var(cube[next]);
93 else
94 vv[next] = (float) 0.0;
95
96 if (cube[i].vol > 1)
97 vv[i] = Var(cube[i]);
98 else
99 vv[i] = (float) 0.0;
100 } else {
101 vv[next] = (float) 0.0;
102 i--;
103 }
104 next = 0;
105 temp = vv[0];
106 for (k = 1; k <= i; k++)
107 if (vv[k] > temp) {
108 temp = vv[k];
109 next = k;
110 }
111 if (temp <= 0.0) {
112 K = i + 1;
113 System.out.println("Only got " + K + " boxes.");
114 break;
115 }
116 }
117 System.out.println("Partition done...");
118
119
120 for (k = 0; k < K; k++) {
121 Mark(cube[k], k, tag);
122
123 weight = Vol(cube[k], wt);
124 if (weight != 0) {
125 lutr[k] = (short) (Vol(cube[k], mr) / weight);
126 lutg[k] = (short) (Vol(cube[k], mg) / weight);
127 lutb[k] = (short) (Vol(cube[k], mb) / weight);
128 } else {
129 System.out.println("Bogus box " + k);
130 lutr[k] = lutg[k] = lutb[k] = 0;
131 }
132 }
133
134
135 for (y = 0; y < h; y++)
136 for (x = 0; x < w; x++) {
137 int pr = qadd[x][y] / 1089;
138 int pt = qadd[x][y] % 1089;
139 int pg = pt / 33;
140 int pb = pt % 33;
141
142 qadd[x][y] = tag[pr][pg][pb];
143 }
144
145
146 for (y = 0; y < h; y++)
147 for (x = 0; x < w; x++) {
148 r[x][y] = lutr[qadd[x][y]];
149 g[x][y] = lutg[qadd[x][y]];
150 b[x][y] = lutb[qadd[x][y]];
151 }
152
153 }
154
155 public void Hist3d(int vwt[][][],
156 int vmr[][][],
157 int vmg[][][],
158 int vmb[][][],
159 float vm2[][][]) {
160 short red, green, blue;
161 int i, inr, ing, inb;
162 int table[] = new int[MAXCOLOR];
163 short tr, tg, tb;
164
165 for (i = 0; i < MAXCOLOR; i++)
166 table[i] = i * i;
167
168 qadd = new int[width][height];
169
170 for (int y = 0; y < height; y++)
171 for (int x = 0; x < width; x++) {
172 tr = ir[x][y];
173 tg = ig[x][y];
174 tb = ib[x][y];
175 inr = (tr >> 3) + 1;
176 ing = (tg >> 3) + 1;
177 inb = (tb >> 3) + 1;
178 qadd[x][y] = ((inr << 10) +
179 (inr << 6) +
180 inr +
181 (ing << 5) +
182 ing +
183 inb);
184 vwt[inr][ing][inb]++;
185 vmr[inr][ing][inb] += tr;
186 vmg[inr][ing][inb] += tg;
187 vmb[inr][ing][inb] += tb;
188 vm2[inr][ing][inb] +=
189 (float) (table[tr] + table[tg] + table[tb]);
190 }
191 }
192
193 public void M3d(int vwt[][][],
194 int vmr[][][],
195 int vmg[][][],
196 int vmb[][][],
197 float vm2[][][]) {
198 short i, tr, tg, tb;
199 int line, liner, lineg, lineb;
200 int area[] = new int[33];
201 int arear[] = new int[33];
202 int areag[] = new int[33];
203 int areab[] = new int[33];
204 float line2;
205 float area2[] = new float[33];
206
207 for (tr = 1; tr <= 32; tr++) {
208 for (i = 0; i <= 32; i++)
209 area2[i] = area[i] = arear[i] = areag[i] = areab[i] = 0;
210 for (tg = 1; tg <= 32; tg++) {
211 line2 = line = liner = lineg = lineb = 0;
212 for (tb = 1; tb <= 32; tb++) {
213 line += vwt[tr][tg][tb];
214 liner += vmr[tr][tg][tb];
215 lineg += vmg[tr][tg][tb];
216 lineb += vmb[tr][tg][tb];
217 line2 += vm2[tr][tg][tb];
218
219 area[tb] += line;
220 arear[tb] += liner;
221 areag[tb] += lineg;
222 areab[tb] += lineb;
223 area2[tb] += line2;
224
225 vwt[tr][tg][tb] = vwt[tr - 1][tg][tb] + area[tb];
226 vmr[tr][tg][tb] = vmr[tr - 1][tg][tb] + arear[tb];
227 vmg[tr][tg][tb] = vmg[tr - 1][tg][tb] + areag[tb];
228 vmb[tr][tg][tb] = vmb[tr - 1][tg][tb] + areab[tb];
229 vm2[tr][tg][tb] = vm2[tr - 1][tg][tb] + area2[tb];
230
231 }
232 }
233 }
234
235 }
236
237 public long Vol(Box cube, int mmt[][][]) {
238 return (mmt[cube.r1][cube.g1][cube.b1]
239 -
240 mmt[cube.r1][cube.g1][cube.b0]
241 - mmt[cube.r1][cube.g0][cube.b1]
242 + mmt[cube.r1][cube.g0][cube.b0]
243 - mmt[cube.r0][cube.g1][cube.b1]
244 +
245 mmt[cube.r0][cube.g1][cube.b0]
246 + mmt[cube.r0][cube.g0][cube.b1]
247 - mmt[cube.r0][cube.g0][cube.b0]);
248 }
249
250 public long Bottom(Box cube, int dir, int mmt[][][]) {
251 switch (dir) {
252 case RED:
253 return (-mmt[cube.r0][cube.g1][cube.b1]
254 +
255 mmt[cube.r0][cube.g1][cube.b0]
256 + mmt[cube.r0][cube.g0][cube.b1]
257 - mmt[cube.r0][cube.g0][cube.b0]);
258 case GREEN:
259 return (-mmt[cube.r1][cube.g0][cube.b1]
260 +
261 mmt[cube.r1][cube.g0][cube.b0]
262 + mmt[cube.r0][cube.g0][cube.b1]
263 - mmt[cube.r0][cube.g0][cube.b0]);
264 case BLUE:
265 return (-mmt[cube.r1][cube.g1][cube.b0]
266 +
267 mmt[cube.r1][cube.g0][cube.b0]
268 + mmt[cube.r0][cube.g1][cube.b0]
269 - mmt[cube.r0][cube.g0][cube.b0]);
270 }
271
272 return (long) 1;
273 }
274
275 public long Top(Box cube, int dir, int pos,
276 int mmt[][][]) {
277 switch (dir) {
278 case RED:
279 return (mmt[pos][cube.g1][cube.b1]
280 -
281 mmt[pos][cube.g1][cube.b0]
282 - mmt[pos][cube.g0][cube.b1]
283 + mmt[pos][cube.g0][cube.b0]);
284 case GREEN:
285 return (mmt[cube.r1][pos][cube.b1]
286 -
287 mmt[cube.r1][pos][cube.b0]
288 - mmt[cube.r0][pos][cube.b1]
289 + mmt[cube.r0][pos][cube.b0]);
290 case BLUE:
291 return (mmt[cube.r1][cube.g1][pos]
292 -
293 mmt[cube.r1][cube.g0][pos]
294 - mmt[cube.r0][cube.g1][pos]
295 + mmt[cube.r0][cube.g0][pos]);
296 }
297
298 return (long) 1;
299 }
300
301 public float Var(Box cube) {
302 float dr, dg, db, xx;
303
304 dr = Vol(cube, mr);
305 dg = Vol(cube, mg);
306 db = Vol(cube, mb);
307 xx = m2[cube.r1][cube.g1][cube.b1]
308 -
309 m2[cube.r1][cube.g1][cube.b0]
310 - m2[cube.r1][cube.g0][cube.b1]
311 + m2[cube.r1][cube.g0][cube.b0]
312 - m2[cube.r0][cube.g1][cube.b1]
313 +
314 m2[cube.r0][cube.g1][cube.b0]
315 + m2[cube.r0][cube.g0][cube.b1]
316 - m2[cube.r0][cube.g0][cube.b0];
317
318 return (xx -
319 (dr * dr + dg * dg + db * db) / (float) Vol(cube, wt));
320 }
321
322 public float Maximize(Box cube, int dir, int first, int last,
323 int cut[],
324 long wholer,
325 long wholeg,
326 long wholeb,
327 long wholew) {
328 long halfr, halfg, halfb, halfw;
329 long baser, baseg, baseb, basew;
330 int i;
331 float temp, max;
332
333 baser = Bottom(cube, dir, mr);
334 baseg = Bottom(cube, dir, mg);
335 baseb = Bottom(cube, dir, mb);
336 basew = Bottom(cube, dir, wt);
337 max = (float) 0.0;
338 cut[0] = -1;
339
340 for (i = first; i < last; i++) {
341 halfr = baser + Top(cube, dir, i, mr);
342 halfg = baseg + Top(cube, dir, i, mg);
343 halfb = baseb + Top(cube, dir, i, mb);
344 halfw = basew + Top(cube, dir, i, wt);
345
346 if (halfw == 0)
347 continue;
348 else
349 temp = ((float) halfr * halfr +
350 (float) halfg * halfg +
351 (float) halfb * halfb) / halfw;
352
353 halfr = wholer - halfr;
354 halfg = wholeg - halfg;
355 halfb = wholeb - halfb;
356 halfw = wholew - halfw;
357
358 if (halfw == 0)
359 continue;
360 else
361 temp += ((float) halfr * halfr +
362 (float) halfg * halfg +
363 (float) halfb * halfb) / halfw;
364
365 if (temp > max) {
366 max = temp;
367 cut[0] = i;
368 }
369 }
370
371 return (max);
372 }
373
374 public int Cut(Box set1, Box set2) {
375 int dir;
376 int cutr[] = new int[1];
377 int cutg[] = new int[1];
378 int cutb[] = new int[1];
379 float maxr, maxg, maxb;
380 long wholer, wholeg, wholeb, wholew;
381
382 wholer = Vol(set1, mr);
383 wholeg = Vol(set1, mg);
384 wholeb = Vol(set1, mb);
385 wholew = Vol(set1, wt);
386
387
388 maxr = Maximize(set1, RED, set1.r0 + 1, set1.r1, cutr,
389 wholer, wholeg, wholeb, wholew);
390 maxg = Maximize(set1, GREEN, set1.g0 + 1, set1.g1, cutg,
391 wholer, wholeg, wholeb, wholew);
392 maxb = Maximize(set1, BLUE, set1.b0 + 1, set1.b1, cutb,
393 wholer, wholeg, wholeb, wholew);
394 if ((maxr >= maxg) && (maxr >= maxg)) {
395 dir = RED;
396 if (cutr[0] < 0)
397 return 0;
398 } else if ((maxg >= maxr) && (maxg >= maxb)) {
399 dir = GREEN;
400 } else
401 dir = BLUE;
402
403
404 set2.r1 = set1.r1;
405 set2.g1 = set1.g1;
406 set2.b1 = set1.b1;
407
408 switch (dir) {
409 case RED:
410 set2.r0 = set1.r1 = cutr[0];
411 set2.g0 = set1.g0;
412 set2.b0 = set1.b0;
413 break;
414 case GREEN:
415 set2.r0 = set1.r0;
416 set2.g0 = set1.g1 = cutg[0];
417 set2.b0 = set1.b0;
418 break;
419 case BLUE:
420 set2.r0 = set1.r0;
421 set2.g0 = set1.g0;
422 set2.b0 = set1.b1 = cutb[0];
423 break;
424 }
425
426 set1.vol = (set1.r1 - set1.r0) *
427 (set1.g1 - set1.g0) *
428 (set1.b1 - set1.b0);
429 set2.vol = (set2.r1 - set2.r0) *
430 (set2.g1 - set2.g0) *
431 (set2.b1 - set2.b0);
432 return 1;
433 }
434
435 public void Mark(Box cube, int label, int tag[][][]) {
436 int tr, tg, tb;
437
438 for (tr = cube.r0 + 1; tr <= cube.r1; tr++)
439 for (tg = cube.g0 + 1; tg <= cube.g1; tg++)
440 for (tb = cube.b0 + 1; tb <= cube.b1; tb++)
441 tag[tr][tg][tb] = label;
442 }
443 }
444