1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include <float.h>
6
7 /*!
8 **
9 ** Copyright (c) 20011 by John W. Ratcliff mailto:jratcliffscarab@gmail.com
10 **
11 **
12 ** The MIT license:
13 **
14 ** Permission is hereby granted, free of charge, to any person obtaining a copy
15 ** of this software and associated documentation files (the "Software"), to deal
16 ** in the Software without restriction, including without limitation the rights
17 ** to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
18 ** copies of the Software, and to permit persons to whom the Software is furnished
19 ** to do so, subject to the following conditions:
20 **
21 ** The above copyright notice and this permission notice shall be included in all
22 ** copies or substantial portions of the Software.
23
24 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
25 ** IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 ** FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
27 ** AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
28 ** WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
29 ** CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
30
31 */
32
33 #include "WuQuantizer.h"
34 #include "SparseArray.h"
35
36 ///////////////////////////////////////////////////////////////////////
37 // C Implementation of Wu's Color Quantizer (v. 2)
38 // (see Graphics Gems vol. II, pp. 126-133)
39 //
40 // Author: Xiaolin Wu
41 // Dept. of Computer Science
42 // Univ. of Western Ontario
43 // London, Ontario N6A 5B7
44 // wu@csd.uwo.ca
45 //
46 // Algorithm: Greedy orthogonal bipartition of RGB space for variance
47 // minimization aided by inclusion-exclusion tricks.
48 // For speed no nearest neighbor search is done. Slightly
49 // better performance can be expected by more sophisticated
50 // but more expensive versions.
51 //
52 // The author thanks Tom Lane at Tom_Lane@G.GP.CS.CMU.EDU for much of
53 // additional documentation and a cure to a previous bug.
54 //
55 // Free to distribute, comments and suggestions are appreciated.
56 ///////////////////////////////////////////////////////////////////////
57
58 ///////////////////////////////////////////////////////////////////////
59 // History
60 // -------
61 // July 2000: C++ Implementation of Wu's Color Quantizer
62 // and adaptation for the FreeImage 2 Library
63 // Author: Herv� Drolon (drolon@infonie.fr)
64 // March 2004: Adaptation for the FreeImage 3 library (port to big endian processors)
65 // Author: Herv� Drolon (drolon@infonie.fr)
66 ///////////////////////////////////////////////////////////////////////
67
68
69 ///////////////////////////////////////////////////////////////////////
70
71
72
73 #pragma warning(disable:4100)
74
75 using namespace hacd;
76
77 namespace HACD
78 {
79
80 #define FI_RGBA_RED 2
81 #define FI_RGBA_GREEN 1
82 #define FI_RGBA_BLUE 0
83 #define FI_RGBA_ALPHA 3
84 #define FI_RGBA_RED_MASK 0x00FF0000
85 #define FI_RGBA_GREEN_MASK 0x0000FF00
86 #define FI_RGBA_BLUE_MASK 0x000000FF
87 #define FI_RGBA_ALPHA_MASK 0xFF000000
88 #define FI_RGBA_RED_SHIFT 16
89 #define FI_RGBA_GREEN_SHIFT 8
90 #define FI_RGBA_BLUE_SHIFT 0
91 #define FI_RGBA_ALPHA_SHIFT 24
92
93 ////////////////////////////////////////////////////////////////
94 class Vec3
95 {
96 public:
97
Vec3(HaF32 _x,HaF32 _y,HaF32 _z)98 Vec3(HaF32 _x,HaF32 _y,HaF32 _z)
99 {
100 x = _x;
101 y = _y;
102 z = _z;
103 }
104
Vec3(void)105 Vec3(void)
106 {
107 }
108
Vec3(const hacd::HaF32 * p)109 Vec3(const hacd::HaF32 *p)
110 {
111 x = p[0];
112 y = p[1];
113 z = p[2];
114 }
115
116 HaF32 x;
117 HaF32 y;
118 HaF32 z;
119 };
120
121 typedef hacd::vector< Vec3 > Vec3Vector;
122
123
124 /**
125 Xiaolin Wu color quantization algorithm
126 */
127 class WuColorQuantizer
128 {
129 public:
130 // Constructor - Input parameter: DIB 24-bit to be quantized
131 WuColorQuantizer(void);
132 // Destructor
133 ~WuColorQuantizer();
134 // Quantizer - Return value: quantized 8-bit (color palette) DIB
135 void Quantize(HaI32 PaletteSize,Vec3Vector &outputVertices);
136
137 void addColor(HaF32 x,HaF32 y,HaF32 z);
138
139
140 typedef struct tagBox
141 {
142 HaI32 r0; // min value, exclusive
143 HaI32 r1; // max value, inclusive
144 HaI32 g0;
145 HaI32 g1;
146 HaI32 b0;
147 HaI32 b1;
148 HaI32 vol;
149 } Box;
150
151 private:
152 HaI32 table[256];
153 HaF32 *mSumSquared;
154 HaI32 *mWeight;
155 HaI32 *mSumX;
156 HaI32 *mSumY;
157 HaI32 *mSumZ;
158
159 void M3D(HaI32 *vwt, HaI32 *vmr, HaI32 *vmg, HaI32 *vmb, HaF32 *m2);
160 HaI32 Vol(Box *cube, HaI32 *mmt);
161 HaI32 Bottom(Box *cube, HaU8 dir, HaI32 *mmt);
162 HaI32 Top(Box *cube, HaU8 dir, HaI32 pos, HaI32 *mmt);
163 HaF32 Var(Box *cube);
164 HaF32 Maximize(Box *cube, HaU8 dir, HaI32 first, HaI32 last , HaI32 *cut,
165 HaI32 whole_r, HaI32 whole_g, HaI32 whole_b, HaI32 whole_w);
166 bool Cut(Box *set1, Box *set2);
167 void Mark(Box *cube, HaI32 label, HaU8 *tag);
168
169 };
170
171
172
173 // Size of a 3D array : 33 x 33 x 33
174 #define SIZE_3D 35937
175
176 // 3D array indexation
177 #define INDEX(r, g, b) ((r << 10) + (r << 6) + r + (g << 5) + g + b)
178
179 #define MAXCOLOR 1024
180
181 // Constructor / Destructor
182
WuColorQuantizer(void)183 WuColorQuantizer::WuColorQuantizer(void)
184 {
185 // Allocate 3D arrays
186 mSumSquared = (HaF32*)malloc(SIZE_3D * sizeof(HaF32));
187 mWeight = (HaI32*)malloc(SIZE_3D * sizeof(HaI32));
188 mSumX = (HaI32*)malloc(SIZE_3D * sizeof(HaI32));
189 mSumY = (HaI32*)malloc(SIZE_3D * sizeof(HaI32));
190 mSumZ = (HaI32*)malloc(SIZE_3D * sizeof(HaI32));
191
192 memset(mSumSquared, 0, SIZE_3D * sizeof(HaF32));
193 memset(mWeight, 0, SIZE_3D * sizeof(HaI32));
194 memset(mSumX, 0, SIZE_3D * sizeof(HaI32));
195 memset(mSumY, 0, SIZE_3D * sizeof(HaI32));
196 memset(mSumZ, 0, SIZE_3D * sizeof(HaI32));
197
198 for(HaI32 i = 0; i < 256; i++)
199 table[i] = i * i;
200
201 }
202
~WuColorQuantizer()203 WuColorQuantizer::~WuColorQuantizer()
204 {
205 if(mSumSquared) free(mSumSquared);
206 if(mWeight) free(mWeight);
207 if(mSumX) free(mSumX);
208 if(mSumY) free(mSumY);
209 if(mSumZ) free(mSumZ);
210 }
211
212
213 // Histogram is in elements 1..HISTSIZE along each axis,
214 // element 0 is for base or marginal value
215 // NB: these must start out 0!
216
217 // Build 3-D color histogram of counts, r/g/b, c^2
218
addColor(HaF32 x,HaF32 y,HaF32 z)219 void WuColorQuantizer::addColor(HaF32 x,HaF32 y,HaF32 z)
220 {
221 HaU32 red = (HaU32)(x*128+128);
222 HaU32 green = (HaU32)(y*128+128);
223 HaU32 blue = (HaU32)(z*128+128);
224 HACD_ASSERT( red < 256 );
225 HACD_ASSERT( green < 256 );
226 HACD_ASSERT( blue < 256 );
227
228 HaU32 inr = (red >> 3) + 1;
229 HaU32 ing = (green >> 3) + 1;
230 HaU32 inb = (blue >> 3) + 1;
231 HaU32 ind = INDEX(inr, ing, inb);
232
233 mWeight[ind]++;
234 mSumX[ind] += red;
235 mSumY[ind] += green;
236 mSumZ[ind] += blue;
237 mSumSquared[ind] += table[red] + table[green] + table[blue];
238 }
239
240
241 // At conclusion of the histogram step, we can interpret
242 // mWeight[r][g][b] = sum over voxel of P(c)
243 // mSumX[r][g][b] = sum over voxel of r*P(c) , similarly for mSumY, mSumZ
244 // m2[r][g][b] = sum over voxel of c^2*P(c)
245 // Actually each of these should be divided by 'ImageSize' to give the usual
246 // interpretation of P() as ranging from 0 to 1, but we needn't do that here.
247
248
249 // We now convert histogram into moments so that we can rapidly calculate
250 // the sums of the above quantities over any desired box.
251
252 // Compute cumulative moments
M3D(HaI32 * vwt,HaI32 * vmr,HaI32 * vmg,HaI32 * vmb,HaF32 * m2)253 void WuColorQuantizer::M3D(HaI32 *vwt, HaI32 *vmr, HaI32 *vmg, HaI32 *vmb, HaF32 *m2)
254 {
255 HaU32 ind1, ind2;
256 HaU8 i, r, g, b;
257 HaI32 line, line_r, line_g, line_b;
258 HaI32 area[33], area_r[33], area_g[33], area_b[33];
259 HaF32 line2, area2[33];
260
261 for(r = 1; r <= 32; r++)
262 {
263 for(i = 0; i <= 32; i++)
264 {
265 area2[i] = 0;
266 area[i] = area_r[i] = area_g[i] = area_b[i] = 0;
267 }
268 for(g = 1; g <= 32; g++)
269 {
270 line2 = 0;
271 line = line_r = line_g = line_b = 0;
272 for(b = 1; b <= 32; b++)
273 {
274 ind1 = INDEX(r, g, b); // [r][g][b]
275 line += vwt[ind1];
276 line_r += vmr[ind1];
277 line_g += vmg[ind1];
278 line_b += vmb[ind1];
279 line2 += m2[ind1];
280 area[b] += line;
281 area_r[b] += line_r;
282 area_g[b] += line_g;
283 area_b[b] += line_b;
284 area2[b] += line2;
285 ind2 = ind1 - 1089; // [r-1][g][b]
286 vwt[ind1] = vwt[ind2] + area[b];
287 vmr[ind1] = vmr[ind2] + area_r[b];
288 vmg[ind1] = vmg[ind2] + area_g[b];
289 vmb[ind1] = vmb[ind2] + area_b[b];
290 m2[ind1] = m2[ind2] + area2[b];
291 }
292 }
293 }
294 }
295
296 // Compute sum over a box of any given statistic
297 HaI32
Vol(Box * cube,HaI32 * mmt)298 WuColorQuantizer::Vol( Box *cube, HaI32 *mmt ) {
299 return( mmt[INDEX(cube->r1, cube->g1, cube->b1)]
300 - mmt[INDEX(cube->r1, cube->g1, cube->b0)]
301 - mmt[INDEX(cube->r1, cube->g0, cube->b1)]
302 + mmt[INDEX(cube->r1, cube->g0, cube->b0)]
303 - mmt[INDEX(cube->r0, cube->g1, cube->b1)]
304 + mmt[INDEX(cube->r0, cube->g1, cube->b0)]
305 + mmt[INDEX(cube->r0, cube->g0, cube->b1)]
306 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
307 }
308
309 // The next two routines allow a slightly more efficient calculation
310 // of Vol() for a proposed subbox of a given box. The sum of Top()
311 // and Bottom() is the Vol() of a subbox split in the given direction
312 // and with the specified new upper bound.
313
314
315 // Compute part of Vol(cube, mmt) that doesn't depend on r1, g1, or b1
316 // (depending on dir)
317
318 HaI32
Bottom(Box * cube,HaU8 dir,HaI32 * mmt)319 WuColorQuantizer::Bottom(Box *cube, HaU8 dir, HaI32 *mmt) {
320 switch(dir)
321 {
322 case FI_RGBA_RED:
323 return( - mmt[INDEX(cube->r0, cube->g1, cube->b1)]
324 + mmt[INDEX(cube->r0, cube->g1, cube->b0)]
325 + mmt[INDEX(cube->r0, cube->g0, cube->b1)]
326 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
327 break;
328 case FI_RGBA_GREEN:
329 return( - mmt[INDEX(cube->r1, cube->g0, cube->b1)]
330 + mmt[INDEX(cube->r1, cube->g0, cube->b0)]
331 + mmt[INDEX(cube->r0, cube->g0, cube->b1)]
332 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
333 break;
334 case FI_RGBA_BLUE:
335 return( - mmt[INDEX(cube->r1, cube->g1, cube->b0)]
336 + mmt[INDEX(cube->r1, cube->g0, cube->b0)]
337 + mmt[INDEX(cube->r0, cube->g1, cube->b0)]
338 - mmt[INDEX(cube->r0, cube->g0, cube->b0)] );
339 break;
340 }
341
342 return 0;
343 }
344
345
346 // Compute remainder of Vol(cube, mmt), substituting pos for
347 // r1, g1, or b1 (depending on dir)
348
Top(Box * cube,HaU8 dir,HaI32 pos,HaI32 * mmt)349 HaI32 WuColorQuantizer::Top(Box *cube, HaU8 dir, HaI32 pos, HaI32 *mmt)
350 {
351 switch(dir)
352 {
353 case FI_RGBA_RED:
354 return( mmt[INDEX(pos, cube->g1, cube->b1)]
355 -mmt[INDEX(pos, cube->g1, cube->b0)]
356 -mmt[INDEX(pos, cube->g0, cube->b1)]
357 +mmt[INDEX(pos, cube->g0, cube->b0)] );
358 break;
359 case FI_RGBA_GREEN:
360 return( mmt[INDEX(cube->r1, pos, cube->b1)]
361 -mmt[INDEX(cube->r1, pos, cube->b0)]
362 -mmt[INDEX(cube->r0, pos, cube->b1)]
363 +mmt[INDEX(cube->r0, pos, cube->b0)] );
364 break;
365 case FI_RGBA_BLUE:
366 return( mmt[INDEX(cube->r1, cube->g1, pos)]
367 -mmt[INDEX(cube->r1, cube->g0, pos)]
368 -mmt[INDEX(cube->r0, cube->g1, pos)]
369 +mmt[INDEX(cube->r0, cube->g0, pos)] );
370 break;
371 }
372
373 return 0;
374 }
375
376 // Compute the weighted variance of a box
377 // NB: as with the raw statistics, this is really the variance * ImageSize
378
Var(Box * cube)379 HaF32 WuColorQuantizer::Var(Box *cube)
380 {
381 HaF32 dr = (HaF32) Vol(cube, mSumX);
382 HaF32 dg = (HaF32) Vol(cube, mSumY);
383 HaF32 db = (HaF32) Vol(cube, mSumZ);
384 HaF32 xx = mSumSquared[INDEX(cube->r1, cube->g1, cube->b1)]
385 -mSumSquared[INDEX(cube->r1, cube->g1, cube->b0)]
386 -mSumSquared[INDEX(cube->r1, cube->g0, cube->b1)]
387 +mSumSquared[INDEX(cube->r1, cube->g0, cube->b0)]
388 -mSumSquared[INDEX(cube->r0, cube->g1, cube->b1)]
389 +mSumSquared[INDEX(cube->r0, cube->g1, cube->b0)]
390 +mSumSquared[INDEX(cube->r0, cube->g0, cube->b1)]
391 -mSumSquared[INDEX(cube->r0, cube->g0, cube->b0)];
392
393 return (xx - (dr*dr+dg*dg+db*db)/(HaF32)Vol(cube,mWeight));
394 }
395
396 // We want to minimize the sum of the variances of two subboxes.
397 // The sum(c^2) terms can be ignored since their sum over both subboxes
398 // is the same (the sum for the whole box) no matter where we split.
399 // The remaining terms have a minus sign in the variance formula,
400 // so we drop the minus sign and MAXIMIZE the sum of the two terms.
401
Maximize(Box * cube,HaU8 dir,HaI32 first,HaI32 last,HaI32 * cut,HaI32 whole_r,HaI32 whole_g,HaI32 whole_b,HaI32 whole_w)402 HaF32 WuColorQuantizer::Maximize(Box *cube, HaU8 dir, HaI32 first, HaI32 last , HaI32 *cut, HaI32 whole_r, HaI32 whole_g, HaI32 whole_b, HaI32 whole_w)
403 {
404 HaI32 half_r, half_g, half_b, half_w;
405 HaI32 i;
406 HaF32 temp;
407
408 HaI32 base_r = Bottom(cube, dir, mSumX);
409 HaI32 base_g = Bottom(cube, dir, mSumY);
410 HaI32 base_b = Bottom(cube, dir, mSumZ);
411 HaI32 base_w = Bottom(cube, dir, mWeight);
412
413 HaF32 max = 0.0;
414
415 *cut = -1;
416
417 for (i = first; i < last; i++) {
418 half_r = base_r + Top(cube, dir, i, mSumX);
419 half_g = base_g + Top(cube, dir, i, mSumY);
420 half_b = base_b + Top(cube, dir, i, mSumZ);
421 half_w = base_w + Top(cube, dir, i, mWeight);
422
423 // now half_x is sum over lower half of box, if split at i
424
425 if (half_w == 0) { // subbox could be empty of pixels!
426 continue; // never split into an empty box
427 } else {
428 temp = ((HaF32)half_r*half_r + (HaF32)half_g*half_g + (HaF32)half_b*half_b)/half_w;
429 }
430
431 half_r = whole_r - half_r;
432 half_g = whole_g - half_g;
433 half_b = whole_b - half_b;
434 half_w = whole_w - half_w;
435
436 if (half_w == 0) { // subbox could be empty of pixels!
437 continue; // never split into an empty box
438 } else {
439 temp += ((HaF32)half_r*half_r + (HaF32)half_g*half_g + (HaF32)half_b*half_b)/half_w;
440 }
441
442 if (temp > max) {
443 max=temp;
444 *cut=i;
445 }
446 }
447
448 return max;
449 }
450
451 bool
Cut(Box * set1,Box * set2)452 WuColorQuantizer::Cut(Box *set1, Box *set2) {
453 HaU8 dir;
454 HaI32 cutr, cutg, cutb;
455
456 HaI32 whole_r = Vol(set1, mSumX);
457 HaI32 whole_g = Vol(set1, mSumY);
458 HaI32 whole_b = Vol(set1, mSumZ);
459 HaI32 whole_w = Vol(set1, mWeight);
460
461 HaF32 maxr = Maximize(set1, FI_RGBA_RED, set1->r0+1, set1->r1, &cutr, whole_r, whole_g, whole_b, whole_w);
462 HaF32 maxg = Maximize(set1, FI_RGBA_GREEN, set1->g0+1, set1->g1, &cutg, whole_r, whole_g, whole_b, whole_w);
463 HaF32 maxb = Maximize(set1, FI_RGBA_BLUE, set1->b0+1, set1->b1, &cutb, whole_r, whole_g, whole_b, whole_w);
464
465 if ((maxr >= maxg) && (maxr >= maxb)) {
466 dir = FI_RGBA_RED;
467
468 if (cutr < 0) {
469 return false; // can't split the box
470 }
471 } else if ((maxg >= maxr) && (maxg>=maxb)) {
472 dir = FI_RGBA_GREEN;
473 } else {
474 dir = FI_RGBA_BLUE;
475 }
476
477 set2->r1 = set1->r1;
478 set2->g1 = set1->g1;
479 set2->b1 = set1->b1;
480
481 switch (dir) {
482 case FI_RGBA_RED:
483 set2->r0 = set1->r1 = cutr;
484 set2->g0 = set1->g0;
485 set2->b0 = set1->b0;
486 break;
487
488 case FI_RGBA_GREEN:
489 set2->g0 = set1->g1 = cutg;
490 set2->r0 = set1->r0;
491 set2->b0 = set1->b0;
492 break;
493
494 case FI_RGBA_BLUE:
495 set2->b0 = set1->b1 = cutb;
496 set2->r0 = set1->r0;
497 set2->g0 = set1->g0;
498 break;
499 }
500
501 set1->vol = (set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0);
502 set2->vol = (set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0);
503
504 return true;
505 }
506
507
508 void
Mark(Box * cube,HaI32 label,HaU8 * tag)509 WuColorQuantizer::Mark(Box *cube, HaI32 label, HaU8 *tag) {
510 for (HaI32 r = cube->r0 + 1; r <= cube->r1; r++) {
511 for (HaI32 g = cube->g0 + 1; g <= cube->g1; g++) {
512 for (HaI32 b = cube->b0 + 1; b <= cube->b1; b++) {
513 tag[INDEX(r, g, b)] = (HaU8)label;
514 }
515 }
516 }
517 }
518
519 // Wu Quantization algorithm
Quantize(HaI32 PaletteSize,Vec3Vector & outputVertices)520 void WuColorQuantizer::Quantize(HaI32 PaletteSize,Vec3Vector &outputVertices)
521 {
522 HaU8 *tag = NULL;
523
524 if ( PaletteSize > MAXCOLOR )
525 {
526 PaletteSize = MAXCOLOR;
527 }
528 Box cube[MAXCOLOR];
529 HaI32 next;
530 HaI32 i, weight;
531 HaI32 k;
532 HaF32 vv[MAXCOLOR], temp;
533
534 // Compute moments
535 M3D(mWeight, mSumX, mSumY, mSumZ, mSumSquared);
536
537 cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
538 cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
539 next = 0;
540
541 for (i = 1; i < PaletteSize; i++)
542 {
543 if(Cut(&cube[next], &cube[i]))
544 {
545 // volume test ensures we won't try to cut one-cell box
546 vv[next] = (cube[next].vol > 1) ? Var(&cube[next]) : 0;
547 vv[i] = (cube[i].vol > 1) ? Var(&cube[i]) : 0;
548 }
549 else
550 {
551 vv[next] = 0.0; // don't try to split this box again
552 i--; // didn't create box i
553 }
554
555 next = 0; temp = vv[0];
556
557 for (k = 1; k <= i; k++)
558 {
559 if (vv[k] > temp)
560 {
561 temp = vv[k]; next = k;
562 }
563 }
564
565 if (temp <= 0.0)
566 {
567 PaletteSize = i + 1;
568 // Error: "Only got 'PaletteSize' boxes"
569 break;
570 }
571 }
572
573 // Partition done
574 // the space for array mSumSquared can be freed now
575 free(mSumSquared);
576 mSumSquared = NULL;
577
578 // create an optimized palette
579 tag = (HaU8*) malloc(SIZE_3D * sizeof(HaU8));
580 memset(tag, 0, SIZE_3D * sizeof(HaU8));
581
582 for (k = 0; k < PaletteSize ; k++)
583 {
584 Mark(&cube[k], k, tag);
585 weight = Vol(&cube[k], mWeight);
586
587 if (weight)
588 {
589 HaI32 red = (HaI32)(((HaF32)Vol(&cube[k], mSumX) / (HaF32)weight) + 0.5f);
590 HaI32 green = (HaI32)(((HaF32)Vol(&cube[k], mSumY) / (HaF32)weight) + 0.5f);
591 HaI32 blue = (HaI32)(((HaF32)Vol(&cube[k], mSumZ) / (HaF32)weight) + 0.5f);
592 HACD_ASSERT( red >= 0 && red < 256 );
593 HACD_ASSERT( green >= 0 && green < 256 );
594 HACD_ASSERT( blue >= 0 && blue < 256 );
595 Vec3 v;
596 v.x = (red-128.0f)/128.0f;
597 v.y = (green-128.0f)/128.0f;
598 v.z = (blue-128.0f)/128.0f;
599 outputVertices.push_back(v);
600 }
601 else
602 {
603 }
604 }
605 }
606
607
608 hacd::HaU32 kmeans_cluster3d(const hacd::HaF32 *input, // an array of input 3d data points.
609 hacd::HaU32 inputSize, // the number of input data points.
610 hacd::HaU32 clumpCount, // the number of clumps you wish to product.
611 hacd::HaF32 *outputClusters, // The output array of clumps 3d vectors, should be at least 'clumpCount' in size.
612 hacd::HaU32 *outputIndices, // A set of indices which remaps the input vertices to clumps; should be at least 'inputSize'
613 hacd::HaF32 errorThreshold=0.01f, // The error threshold to converge towards before giving up.
614 hacd::HaF32 collapseDistance=0.01f); // distance so small it is not worth bothering to create a new clump.
615
616 template <class Type> struct Vec3d
617 {
distanceSquaredHACD::Vec3d618 inline Type distanceSquared(const Vec3d &a) const
619 {
620 Type dx = x-a.x;
621 Type dy = y-a.y;
622 Type dz = z-a.z;
623 return dx*dx+dy*dy+dz*dz;
624 }
operator +=HACD::Vec3d625 inline void operator+=(const Vec3d &v)
626 {
627 x+=v.x;
628 y+=v.y;
629 z+=v.z;
630 }
operator *=HACD::Vec3d631 inline void operator*=(const Type v)
632 {
633 x*=v;
634 y*=v;
635 z*=v;
636 }
zeroHACD::Vec3d637 inline void zero(void) { x = 0; y = 0; z = 0; };
638
639 Type x;
640 Type y;
641 Type z;
642 };
643
644 template <class Vec,class Type >
kmeans_cluster(const Vec * input,hacd::HaU32 inputCount,hacd::HaU32 clumpCount,Vec * clusters,hacd::HaU32 * outputIndices,Type threshold,Type collapseDistance)645 hacd::HaU32 kmeans_cluster(const Vec *input,
646 hacd::HaU32 inputCount,
647 hacd::HaU32 clumpCount,
648 Vec *clusters,
649 hacd::HaU32 *outputIndices,
650 Type threshold, // controls how long it works to converge towards a least errors solution.
651 Type collapseDistance) // distance between clumps to consider them to be essentially equal.
652 {
653 hacd::HaU32 convergeCount = 64; // maximum number of iterations attempting to converge to a solution..
654 hacd::HaU32 *counts = (hacd::HaU32 *)HACD_ALLOC(sizeof(hacd::HaU32)*clumpCount);
655 Type error=0;
656 if ( inputCount <= clumpCount ) // if the number of input points is less than our clumping size, just return the input points.
657 {
658 clumpCount = inputCount;
659 for (hacd::HaU32 i=0; i<inputCount; i++)
660 {
661 if ( outputIndices )
662 {
663 outputIndices[i] = i;
664 }
665 clusters[i] = input[i];
666 counts[i] = 1;
667 }
668 }
669 else
670 {
671 Vec *centroids = (Vec *)HACD_ALLOC(sizeof(Vec)*clumpCount);
672
673 // Take a sampling of the input points as initial centroid estimates.
674 for (hacd::HaU32 i=0; i<clumpCount; i++)
675 {
676 hacd::HaU32 index = (i*inputCount)/clumpCount;
677 assert( index >= 0 && index < inputCount );
678 clusters[i] = input[index];
679 }
680
681 // Here is the main convergence loop
682 Type old_error = FLT_MAX; // old and initial error estimates are max Type
683 error = FLT_MAX;
684 do
685 {
686 old_error = error; // preserve the old error
687 // reset the counts and centroids to current cluster location
688 for (hacd::HaU32 i=0; i<clumpCount; i++)
689 {
690 counts[i] = 0;
691 centroids[i].zero();
692 }
693 error = 0;
694 // For each input data point, figure out which cluster it is closest too and add it to that cluster.
695 for (hacd::HaU32 i=0; i<inputCount; i++)
696 {
697 Type min_distance = FLT_MAX;
698 // find the nearest clump to this point.
699 for (hacd::HaU32 j=0; j<clumpCount; j++)
700 {
701 Type distance = input[i].distanceSquared( clusters[j] );
702 if ( distance < min_distance )
703 {
704 min_distance = distance;
705 outputIndices[i] = j; // save which clump this point indexes
706 }
707 }
708 hacd::HaU32 index = outputIndices[i]; // which clump was nearest to this point.
709 centroids[index]+=input[i];
710 counts[index]++; // increment the counter indicating how many points are in this clump.
711 error+=min_distance; // save the error accumulation
712 }
713 // Now, for each clump, compute the mean and store the result.
714 for (hacd::HaU32 i=0; i<clumpCount; i++)
715 {
716 if ( counts[i] ) // if this clump got any points added to it...
717 {
718 Type recip = 1.0f / (Type)counts[i]; // compute the average (center of those points)
719 centroids[i]*=recip; // compute the average center of the points in this clump.
720 clusters[i] = centroids[i]; // store it as the new cluster.
721 }
722 }
723 // decrement the convergence counter and bail if it is taking too long to converge to a solution.
724 convergeCount--;
725 if (convergeCount == 0 )
726 {
727 break;
728 }
729 if ( error < threshold ) // early exit if our first guess is already good enough (if all input points are the same)
730 break;
731 } while ( fabs(error - old_error) > threshold ); // keep going until the error is reduced by this threshold amount.
732
733 HACD_FREE(centroids);
734 }
735
736 // ok..now we prune the clumps if necessary.
737 // The rules are; first, if a clump has no 'counts' then we prune it as it's unused.
738 // The second, is if the centroid of this clump is essentially the same (based on the distance tolerance)
739 // as an existing clump, then it is pruned and all indices which used to point to it, now point to the one
740 // it is closest too.
741 hacd::HaU32 outCount = 0; // number of clumps output after pruning performed.
742 Type d2 = collapseDistance*collapseDistance; // squared collapse distance.
743 for (hacd::HaU32 i=0; i<clumpCount; i++)
744 {
745 if ( counts[i] == 0 ) // if no points ended up in this clump, eliminate it.
746 continue;
747 // see if this clump is too close to any already accepted clump.
748 bool add = true;
749 hacd::HaU32 remapIndex = outCount; // by default this clump will be remapped to its current index.
750 for (hacd::HaU32 j=0; j<outCount; j++)
751 {
752 Type distance = clusters[i].distanceSquared(clusters[j]);
753 if ( distance < d2 )
754 {
755 remapIndex = j;
756 add = false; // we do not add this clump
757 break;
758 }
759 }
760 // If we have fewer output clumps than input clumps so far, then we need to remap the old indices to the new ones.
761 if ( outputIndices )
762 {
763 if ( outCount != i || !add ) // we need to remap indices! everything that was index 'i' now needs to be remapped to 'outCount'
764 {
765 for (hacd::HaU32 j=0; j<inputCount; j++)
766 {
767 if ( outputIndices[j] == i )
768 {
769 outputIndices[j] = remapIndex; //
770 }
771 }
772 }
773 }
774 if ( add )
775 {
776 clusters[outCount] = clusters[i];
777 outCount++;
778 }
779 }
780 HACD_FREE(counts);
781 clumpCount = outCount;
782 return clumpCount;
783 }
784
kmeans_cluster3d(const hacd::HaF32 * input,hacd::HaU32 inputSize,hacd::HaU32 clumpCount,hacd::HaF32 * outputClusters,hacd::HaU32 * outputIndices,hacd::HaF32 errorThreshold,hacd::HaF32 collapseDistance)785 hacd::HaU32 kmeans_cluster3d(const hacd::HaF32 *input, // an array of input 3d data points.
786 hacd::HaU32 inputSize, // the number of input data points.
787 hacd::HaU32 clumpCount, // the number of clumps you wish to produce
788 hacd::HaF32 *outputClusters, // The output array of clumps 3d vectors, should be at least 'clumpCount' in size.
789 hacd::HaU32 *outputIndices, // A set of indices which remaps the input vertices to clumps; should be at least 'inputSize'
790 hacd::HaF32 errorThreshold, // The error threshold to converge towards before giving up.
791 hacd::HaF32 collapseDistance) // distance so small it is not worth bothering to create a new clump.
792 {
793 const Vec3d< hacd::HaF32 > *_input = (const Vec3d<hacd::HaF32> *)input;
794 Vec3d<hacd::HaF32> *_output = (Vec3d<hacd::HaF32> *)outputClusters;
795 return kmeans_cluster< Vec3d<hacd::HaF32>, hacd::HaF32 >(_input,inputSize,clumpCount,_output,outputIndices,errorThreshold,collapseDistance);
796 }
797
798 class MyWuQuantizer : public WuQuantizer, public UANS::UserAllocated
799 {
800 public:
MyWuQuantizer(void)801 MyWuQuantizer(void)
802 {
803 mScale = Vec3(1,1,1);
804 mCenter = Vec3(0,0,0);
805 }
806
807 // use the Wu quantizer with 10 bits of resolution on each axis. Precision down to 0.0009765625
808 // All input data is normalized to a unit cube.
809
wuQuantize3D(hacd::HaU32 vcount,const hacd::HaF32 * vertices,bool denormalizeResults,hacd::HaU32 maxVertices,hacd::HaU32 & outputCount)810 virtual const hacd::HaF32 * wuQuantize3D(hacd::HaU32 vcount,
811 const hacd::HaF32 *vertices,
812 bool denormalizeResults,
813 hacd::HaU32 maxVertices,
814 hacd::HaU32 &outputCount)
815 {
816 const hacd::HaF32 *ret = NULL;
817 outputCount = 0;
818
819 normalizeInput(vcount,vertices);
820
821 WuColorQuantizer wcq;
822
823 for (HaU32 i=0; i<vcount; i++)
824 {
825 const Vec3 &v = mNormalizedInput[i];
826 wcq.addColor(v.x,v.y,v.z);
827 }
828
829 wcq.Quantize(maxVertices,mQuantizedOutput);
830
831 outputCount = (HaU32)mQuantizedOutput.size();
832
833 if ( outputCount > 0 )
834 {
835 if ( denormalizeResults )
836 {
837 for (HaU32 i=0; i<outputCount; i++)
838 {
839 Vec3 &v = mQuantizedOutput[i];
840 v.x = v.x*mScale.x + mCenter.x;
841 v.y = v.x*mScale.y + mCenter.y;
842 v.z = v.x*mScale.z + mCenter.z;
843 mQuantizedOutput.push_back(v);
844 }
845 }
846 ret = &mQuantizedOutput[0].x;
847 }
848
849
850 return ret;
851 }
852
853 // Use the kemans quantizer, similar results, but much slower.
kmeansQuantize3D(hacd::HaU32 vcount,const hacd::HaF32 * vertices,bool denormalizeResults,hacd::HaU32 maxVertices,hacd::HaU32 & outputCount)854 virtual const hacd::HaF32 * kmeansQuantize3D(hacd::HaU32 vcount,
855 const hacd::HaF32 *vertices,
856 bool denormalizeResults,
857 hacd::HaU32 maxVertices,
858 hacd::HaU32 &outputCount)
859 {
860 const hacd::HaF32 *ret = NULL;
861 outputCount = 0;
862 mNormalizedInput.clear();
863 mQuantizedOutput.clear();
864
865 if ( vcount > 0 )
866 {
867 normalizeInput(vcount,vertices);
868 HaF32 *quantizedOutput = (HaF32 *)HACD_ALLOC( sizeof(HaF32)*3*vcount);
869 HaU32 *quantizedIndices = (HaU32 *)HACD_ALLOC( sizeof(HaU32)*vcount );
870 outputCount = kmeans_cluster3d(&mNormalizedInput[0].x, vcount, maxVertices, quantizedOutput, quantizedIndices, 0.01f, 0.0001f );
871 if ( outputCount > 0 )
872 {
873 if ( denormalizeResults )
874 {
875 for (HaU32 i=0; i<outputCount; i++)
876 {
877 Vec3 v( &quantizedOutput[i*3] );
878 v.x = v.x*mScale.x + mCenter.x;
879 v.y = v.x*mScale.y + mCenter.y;
880 v.z = v.x*mScale.z + mCenter.z;
881 mQuantizedOutput.push_back(v);
882 }
883 }
884 else
885 {
886 for (HaU32 i=0; i<outputCount; i++)
887 {
888 Vec3 v( &quantizedOutput[i*3] );
889 mQuantizedOutput.push_back(v);
890 }
891 }
892 ret = &mQuantizedOutput[0].x;
893 }
894 HACD_FREE(quantizedOutput);
895 HACD_FREE(quantizedIndices);
896 }
897
898 return ret;
899 }
900
release(void)901 virtual void release(void)
902 {
903 delete this;
904 }
905
getDenormalizeScale(void) const906 virtual const hacd::HaF32 * getDenormalizeScale(void) const
907 {
908 return &mScale.x;
909 }
910
getDenormalizeCenter(void) const911 virtual const hacd::HaF32 * getDenormalizeCenter(void) const
912 {
913 return &mCenter.x;
914 }
915
916
917
918 private:
919
normalizeInput(HaU32 vcount,const HaF32 * vertices)920 void normalizeInput(HaU32 vcount,const HaF32 *vertices)
921 {
922 mNormalizedInput.clear();
923 mQuantizedOutput.clear();
924 Vec3 bmin(vertices);
925 Vec3 bmax(vertices);
926 for (HaU32 i=1; i<vcount; i++)
927 {
928 Vec3 v(&vertices[i*3]);
929
930 if ( v.x < bmin.x )
931 {
932 bmin.x = v.x;
933 }
934 else if ( v.x > bmax.x )
935 {
936 bmax.x = v.x;
937 }
938
939 if ( v.y < bmin.y )
940 {
941 bmin.y = v.y;
942 }
943 else if ( v.y > bmax.y )
944 {
945 bmax.y = v.y;
946 }
947
948 if ( v.z < bmin.z )
949 {
950 bmin.z = v.z;
951 }
952 else if ( v.z > bmax.z )
953 {
954 bmax.z = v.z;
955 }
956 }
957
958 mCenter.x = (bmin.x+bmax.x)*0.5f;
959 mCenter.y = (bmin.y+bmax.y)*0.5f;
960 mCenter.z = (bmin.z+bmax.z)*0.5f;
961
962 HaF32 dx = bmax.x-bmin.x;
963 HaF32 dy = bmax.y-bmin.y;
964 HaF32 dz = bmax.z-bmin.z;
965
966 if ( dx == 0 )
967 {
968 mScale.x = 1;
969 }
970 else
971 {
972 dx = dx*1.001f;
973 mScale.x = dx*0.5f;
974 }
975 if ( dy == 0 )
976 {
977 mScale.y = 1;
978 }
979 else
980 {
981 dy = dy*1.001f;
982 mScale.y = dy*0.5f;
983 }
984 if ( dz == 0 )
985 {
986 mScale.z = 1;
987 }
988 else
989 {
990 dz = dz*1.001f;
991 mScale.z = dz*0.5f;
992 }
993
994 Vec3 recip;
995 recip.x = 1.0f / mScale.x;
996 recip.y = 1.0f / mScale.y;
997 recip.z = 1.0f / mScale.z;
998
999 for (HaU32 i=0; i<vcount; i++)
1000 {
1001 Vec3 v(&vertices[i*3]);
1002
1003 v.x = (v.x-mCenter.x)*recip.x;
1004 v.y = (v.y-mCenter.y)*recip.y;
1005 v.z = (v.z-mCenter.z)*recip.z;
1006
1007 HACD_ASSERT( v.x >= -1 && v.x <= 1 );
1008 HACD_ASSERT( v.y >= -1 && v.y <= 1 );
1009 HACD_ASSERT( v.z >= -1 && v.z <= 1 );
1010 mNormalizedInput.push_back(v);
1011 }
1012 }
1013
~MyWuQuantizer(void)1014 virtual ~MyWuQuantizer(void)
1015 {
1016
1017 }
1018
1019 Vec3 mScale;
1020 Vec3 mCenter;
1021 Vec3Vector mNormalizedInput;
1022 Vec3Vector mQuantizedOutput;
1023 };
1024
createWuQuantizer(void)1025 WuQuantizer * createWuQuantizer(void)
1026 {
1027 MyWuQuantizer *m = HACD_NEW(MyWuQuantizer);
1028 return static_cast< WuQuantizer *>(m);
1029 }
1030
1031
1032 } // end of HACD namespace
1033