← Back to branch summary

~ubuntu-branches/ubuntu/hoary/devil/hoary

« back to all changes in this revision

Viewing changes to src-IL/src/il_quantizer.c

  • Committer: Bazaar Package Importer
  • Author(s): Marcelo E. Magallon
  • Date: 2005-01-03 19:57:42 UTC
  • Revision ID: james.westby@ubuntu.com-20050103195742-4ipkplcwygu3irv0
Tags: upstream-1.6.7
Import upstream version 1.6.7

Show diffs side-by-side

added added

removed removed

Lines of Context:
src-IL/src/il_quantizer.c
 
1
//-----------------------------------------------------------------------
 
2
//         Color Quantization Demo
 
3
//
 
4
// Author: Roman Podobedov
 
5
// Email: romka@ut.ee
 
6
// Romka Graphics: www.ut.ee/~romka
 
7
// 
 
8
// Also in this file implemented Wu's Color Quantizer algorithm (v. 2)
 
9
// For details see Graphics Gems vol. II, pp. 126-133
 
10
//
 
11
// Wu's Color Quantizer Algorithm:
 
12
// Author: Xiaolin Wu
 
13
//                      Dept. of Computer Science
 
14
//                      Univ. of Western Ontario
 
15
//                      London, Ontario N6A 5B7
 
16
//                      wu@csd.uwo.ca
 
17
//                      http://www.csd.uwo.ca/faculty/wu/
 
18
//-----------------------------------------------------------------------
 
19
 
 
20
//-----------------------------------------------------------------------------
 
21
//
 
22
// ImageLib Sources
 
23
// by Denton Woods
 
24
// Last modified: 02/02/2002 <--Y2K Compliant! =]
 
25
//
 
26
// Filename: src-IL/src/il_quantizer.c
 
27
//
 
28
// Description:  Heavily modified by Denton Woods.
 
29
//
 
30
// 20040223 XIX : Modified so it works better with color requests < 256
 
31
// pallete always has memory space for 256 entries
 
32
// used so we can quant down to 255 colors then add a transparent color in there.
 
33
//
 
34
//-----------------------------------------------------------------------------
 
35
 
 
36
#include "il_internal.h"
 
37
 
 
38
#define MAXCOLOR        256
 
39
#define RED                     2
 
40
#define GREEN           1
 
41
#define BLUE            0
 
42
 
 
43
typedef struct Box
 
44
{
 
45
    ILint r0;  // min value, exclusive
 
46
    ILint r1;  // max value, inclusive
 
47
    ILint g0;  
 
48
    ILint g1;  
 
49
    ILint b0;  
 
50
    ILint b1;
 
51
    ILint vol;
 
52
} Box;
 
53
 
 
54
/* Histogram is in elements 1..HISTSIZE along each axis,
 
55
 * element 0 is for base or marginal value
 
56
 * NB: these must start out 0!
 
57
 */
 
58
 
 
59
ILfloat         gm2[33][33][33];
 
60
ILint           wt[33][33][33], mr[33][33][33], mg[33][33][33], mb[33][33][33];
 
61
ILuint          size; //image size
 
62
ILint           K;    //colour look-up table size
 
63
ILushort        *Qadd;
 
64
 
 
65
 
 
66
ILint   WindW, WindH, WindD;
 
67
ILint   i;
 
68
ILubyte *buffer;
 
69
static ILint    Width, Height, Depth, Comp;
 
70
/*ILint TotalColors;
 
71
ILint   a, b;
 
72
ILubyte *buf1, *buf2;*/
 
73
 
 
74
ILuint n2(ILint s)
 
75
{
 
76
        ILint   i;
 
77
        ILint   res;
 
78
 
 
79
        res = 1;
 
80
        for (i = 0; i < s; i++) {
 
81
                res = res*2;
 
82
        }
 
83
 
 
84
        return res;
 
85
}
 
86
 
 
87
 
 
88
// Build 3-D color histogram of counts, r/g/b, c^2
 
89
ILboolean Hist3d(ILubyte *Ir, ILubyte *Ig, ILubyte *Ib, ILint *vwt, ILint *vmr, ILint *vmg, ILint *vmb, ILfloat *m2)
 
90
{
 
91
        ILint   ind, r, g, b;
 
92
        ILint   inr, ing, inb, table[2560];
 
93
        ILuint  i;
 
94
                
 
95
        for (i = 0; i < 256; i++) {
 
96
                table[i] = i * i;
 
97
        }
 
98
        Qadd = (ILushort*)ialloc(sizeof(ILushort) * size);
 
99
        if (Qadd == NULL) {
 
100
                return IL_FALSE;
 
101
        }
 
102
        
 
103
        imemclear(Qadd, sizeof(ILushort) * size);
 
104
        
 
105
        for (i = 0; i < size; i++) {
 
106
            r = Ir[i]; g = Ig[i]; b = Ib[i];
 
107
            inr = (r>>3) + 1; 
 
108
            ing = (g>>3) + 1; 
 
109
            inb = (b>>3) + 1; 
 
110
            Qadd[i] = ind = (inr<<10)+(inr<<6)+inr+(ing<<5)+ing+inb;
 
111
            //[inr][ing][inb]
 
112
            vwt[ind]++;
 
113
            vmr[ind] += r;
 
114
            vmg[ind] += g;
 
115
            vmb[ind] += b;
 
116
                m2[ind] += (ILfloat)(table[r]+table[g]+table[b]);
 
117
        }
 
118
        return IL_TRUE;
 
119
}
 
120
 
 
121
/* At conclusion of the histogram step, we can interpret
 
122
 *   wt[r][g][b] = sum over voxel of P(c)
 
123
 *   mr[r][g][b] = sum over voxel of r*P(c)  ,  similarly for mg, mb
 
124
 *   m2[r][g][b] = sum over voxel of c^2*P(c)
 
125
 * Actually each of these should be divided by 'size' to give the usual
 
126
 * interpretation of P() as ranging from 0 to 1, but we needn't do that here.
 
127
 */
 
128
 
 
129
/* We now convert histogram into moments so that we can rapidly calculate
 
130
 * the sums of the above quantities over any desired Box.
 
131
 */
 
132
 
 
133
 
 
134
// Compute cumulative moments
 
135
ILvoid M3d(ILint *vwt, ILint *vmr, ILint *vmg, ILint *vmb, ILfloat *m2)
 
136
{
 
137
        ILushort        ind1, ind2;
 
138
        ILubyte         i, r, g, b;
 
139
        ILint           line, line_r, line_g, line_b, area[33], area_r[33], area_g[33], area_b[33];
 
140
        ILfloat         line2, area2[33];
 
141
 
 
142
        for (r = 1; r <= 32; r++) {
 
143
                for (i = 0; i <= 32; i++) {
 
144
                        area2[i] = 0.0f;
 
145
                        area[i]=area_r[i]=area_g[i]=area_b[i]=0;
 
146
                }
 
147
                for (g = 1; g <= 32; g++) {
 
148
                        line2 = 0.0f;
 
149
                        line = line_r = line_g = line_b = 0;
 
150
                        for (b = 1; b <= 32; b++) {
 
151
                                ind1 = (r<<10) + (r<<6) + r + (g<<5) + g + b; // [r][g][b]
 
152
                                line += vwt[ind1];
 
153
                                line_r += vmr[ind1]; 
 
154
                                line_g += vmg[ind1]; 
 
155
                                line_b += vmb[ind1];
 
156
                                line2 += m2[ind1];
 
157
                                area[b] += line;
 
158
                                area_r[b] += line_r;
 
159
                                area_g[b] += line_g;
 
160
                                area_b[b] += line_b;
 
161
                                area2[b] += line2;
 
162
                                ind2 = ind1 - 1089; // [r-1][g][b]
 
163
                                vwt[ind1] = vwt[ind2] + area[b];
 
164
                                vmr[ind1] = vmr[ind2] + area_r[b];
 
165
                                vmg[ind1] = vmg[ind2] + area_g[b];
 
166
                                vmb[ind1] = vmb[ind2] + area_b[b];
 
167
                                m2[ind1] = m2[ind2] + area2[b];
 
168
                        }
 
169
                }
 
170
    }
 
171
 
 
172
        return;
 
173
}
 
174
 
 
175
 
 
176
// Compute sum over a Box of any given statistic
 
177
ILint Vol(Box *cube, ILint mmt[33][33][33]) 
 
178
{
 
179
    return( mmt[cube->r1][cube->g1][cube->b1] 
 
180
           -mmt[cube->r1][cube->g1][cube->b0]
 
181
           -mmt[cube->r1][cube->g0][cube->b1]
 
182
           +mmt[cube->r1][cube->g0][cube->b0]
 
183
           -mmt[cube->r0][cube->g1][cube->b1]
 
184
           +mmt[cube->r0][cube->g1][cube->b0]
 
185
           +mmt[cube->r0][cube->g0][cube->b1]
 
186
           -mmt[cube->r0][cube->g0][cube->b0] );
 
187
}
 
188
 
 
189
/* The next two routines allow a slightly more efficient calculation
 
190
 * of Vol() for a proposed subBox of a given Box.  The sum of Top()
 
191
 * and Bottom() is the Vol() of a subBox split in the given direction
 
192
 * and with the specified new upper bound.
 
193
 */
 
194
 
 
195
// Compute part of Vol(cube, mmt) that doesn't depend on r1, g1, or b1
 
196
//      (depending on dir)
 
197
ILint Bottom(Box *cube, ILubyte dir, ILint mmt[33][33][33])
 
198
{
 
199
    switch(dir)
 
200
    {
 
201
                case RED:
 
202
                        return( -mmt[cube->r0][cube->g1][cube->b1]
 
203
                                +mmt[cube->r0][cube->g1][cube->b0]
 
204
                                +mmt[cube->r0][cube->g0][cube->b1]
 
205
                                -mmt[cube->r0][cube->g0][cube->b0] );
 
206
                        break;
 
207
                case GREEN:
 
208
                        return( -mmt[cube->r1][cube->g0][cube->b1]
 
209
                                +mmt[cube->r1][cube->g0][cube->b0]
 
210
                                +mmt[cube->r0][cube->g0][cube->b1]
 
211
                                -mmt[cube->r0][cube->g0][cube->b0] );
 
212
                        break;
 
213
                case BLUE:
 
214
                        return( -mmt[cube->r1][cube->g1][cube->b0]
 
215
                                +mmt[cube->r1][cube->g0][cube->b0]
 
216
                                +mmt[cube->r0][cube->g1][cube->b0]
 
217
                                -mmt[cube->r0][cube->g0][cube->b0] );
 
218
                        break;
 
219
    }
 
220
    return 0;
 
221
}
 
222
 
 
223
 
 
224
// Compute remainder of Vol(cube, mmt), substituting pos for
 
225
//      r1, g1, or b1 (depending on dir)
 
226
ILint Top(Box *cube, ILubyte dir, ILint pos, ILint mmt[33][33][33])
 
227
{
 
228
    switch (dir)
 
229
    {
 
230
                case RED:
 
231
                        return( mmt[pos][cube->g1][cube->b1] 
 
232
                           -mmt[pos][cube->g1][cube->b0]
 
233
                           -mmt[pos][cube->g0][cube->b1]
 
234
                           +mmt[pos][cube->g0][cube->b0] );
 
235
                        break;
 
236
                case GREEN:
 
237
                        return( mmt[cube->r1][pos][cube->b1] 
 
238
                           -mmt[cube->r1][pos][cube->b0]
 
239
                           -mmt[cube->r0][pos][cube->b1]
 
240
                           +mmt[cube->r0][pos][cube->b0] );
 
241
                        break;
 
242
                case BLUE:
 
243
                        return( mmt[cube->r1][cube->g1][pos]
 
244
                           -mmt[cube->r1][cube->g0][pos]
 
245
                           -mmt[cube->r0][cube->g1][pos]
 
246
                           +mmt[cube->r0][cube->g0][pos] );
 
247
                        break;
 
248
    }
 
249
 
 
250
    return 0;
 
251
}
 
252
 
 
253
 
 
254
// Compute the weighted variance of a Box
 
255
//      NB: as with the raw statistics, this is really the variance * size
 
256
ILfloat Var(Box *cube)
 
257
{
 
258
        ILfloat dr, dg, db, xx;
 
259
 
 
260
        dr = (ILfloat)Vol(cube, mr); 
 
261
        dg = (ILfloat)Vol(cube, mg); 
 
262
        db = (ILfloat)Vol(cube, mb);
 
263
        xx = gm2[cube->r1][cube->g1][cube->b1] 
 
264
                -gm2[cube->r1][cube->g1][cube->b0]
 
265
                -gm2[cube->r1][cube->g0][cube->b1]
 
266
                +gm2[cube->r1][cube->g0][cube->b0]
 
267
                -gm2[cube->r0][cube->g1][cube->b1]
 
268
                +gm2[cube->r0][cube->g1][cube->b0]
 
269
                +gm2[cube->r0][cube->g0][cube->b1]
 
270
                -gm2[cube->r0][cube->g0][cube->b0];
 
271
 
 
272
        return xx - (dr*dr+dg*dg+db*db) / (ILfloat)Vol(cube, wt);
 
273
}
 
274
 
 
275
/* We want to minimize the sum of the variances of two subBoxes.
 
276
 * The sum(c^2) terms can be ignored since their sum over both subBoxes
 
277
 * is the same (the sum for the whole Box) no matter where we split.
 
278
 * The remaining terms have a minus sign in the variance formula,
 
279
 * so we drop the minus sign and MAXIMIZE the sum of the two terms.
 
280
 */
 
281
 
 
282
ILfloat Maximize(Box *cube, ILubyte dir, ILint first, ILint last, ILint *cut,
 
283
                                 ILint whole_r, ILint whole_g, ILint whole_b, ILint whole_w)
 
284
{
 
285
        ILint   half_r, half_g, half_b, half_w;
 
286
        ILint   base_r, base_g, base_b, base_w;
 
287
        ILint   i;
 
288
        ILfloat temp, max;
 
289
 
 
290
        base_r = Bottom(cube, dir, mr);
 
291
        base_g = Bottom(cube, dir, mg);
 
292
        base_b = Bottom(cube, dir, mb);
 
293
        base_w = Bottom(cube, dir, wt);
 
294
        max = 0.0;
 
295
        *cut = -1;
 
296
 
 
297
        for (i = first; i < last; ++i) {
 
298
                half_r = base_r + Top(cube, dir, i, mr);
 
299
                half_g = base_g + Top(cube, dir, i, mg);
 
300
                half_b = base_b + Top(cube, dir, i, mb);
 
301
                half_w = base_w + Top(cube, dir, i, wt);
 
302
                // Now half_x is sum over lower half of Box, if split at i 
 
303
                if (half_w == 0) {  // subBox could be empty of pixels!
 
304
                        continue;       // never split into an empty Box
 
305
                }
 
306
                else {
 
307
                        temp = ((ILfloat)half_r*half_r + (ILfloat)half_g * half_g +
 
308
                                        (ILfloat)half_b*half_b) / half_w;
 
309
                }
 
310
 
 
311
                half_r = whole_r - half_r;
 
312
                half_g = whole_g - half_g;
 
313
                half_b = whole_b - half_b;
 
314
                half_w = whole_w - half_w;
 
315
                if (half_w == 0) {  // subBox could be empty of pixels!
 
316
                        continue;       // never split into an empty Box
 
317
                }
 
318
                else {
 
319
                        temp += ((ILfloat)half_r*half_r + (ILfloat)half_g * half_g +
 
320
                                        (ILfloat)half_b*half_b) / half_w;
 
321
                }
 
322
 
 
323
                if (temp > max) {
 
324
                        max = temp;
 
325
                        *cut = i;
 
326
                }
 
327
        }
 
328
 
 
329
        return max;
 
330
}
 
331
 
 
332
 
 
333
ILint Cut(Box *set1, Box *set2)
 
334
{
 
335
        ILubyte dir;
 
336
        ILint cutr, cutg, cutb;
 
337
        ILfloat maxr, maxg, maxb;
 
338
        ILint whole_r, whole_g, whole_b, whole_w;
 
339
 
 
340
        whole_r = Vol(set1, mr);
 
341
        whole_g = Vol(set1, mg);
 
342
        whole_b = Vol(set1, mb);
 
343
        whole_w = Vol(set1, wt);
 
344
 
 
345
        maxr = Maximize(set1, RED, set1->r0+1, set1->r1, &cutr, whole_r, whole_g, whole_b, whole_w);
 
346
        maxg = Maximize(set1, GREEN, set1->g0+1, set1->g1, &cutg, whole_r, whole_g, whole_b, whole_w);
 
347
        maxb = Maximize(set1, BLUE, set1->b0+1, set1->b1, &cutb, whole_r, whole_g, whole_b, whole_w);
 
348
 
 
349
        if ((maxr >= maxg) && (maxr >= maxb)) {
 
350
                dir = RED;
 
351
                if (cutr < 0)
 
352
                        return 0; // can't split the Box 
 
353
        }
 
354
        else if ((maxg >= maxr) && (maxg >= maxb))
 
355
                dir = GREEN;
 
356
        else
 
357
                dir = BLUE;
 
358
 
 
359
        set2->r1 = set1->r1;
 
360
        set2->g1 = set1->g1;
 
361
        set2->b1 = set1->b1;
 
362
 
 
363
        switch (dir)
 
364
        {
 
365
                case RED:
 
366
                        set2->r0 = set1->r1 = cutr;
 
367
                        set2->g0 = set1->g0;
 
368
                        set2->b0 = set1->b0;
 
369
                        break;
 
370
                case GREEN:
 
371
                        set2->g0 = set1->g1 = cutg;
 
372
                        set2->r0 = set1->r0;
 
373
                        set2->b0 = set1->b0;
 
374
                        break;
 
375
                case BLUE:
 
376
                        set2->b0 = set1->b1 = cutb;
 
377
                        set2->r0 = set1->r0;
 
378
                        set2->g0 = set1->g0;
 
379
                        break;
 
380
        }
 
381
 
 
382
        set1->vol = (set1->r1-set1->r0) * (set1->g1-set1->g0) * (set1->b1-set1->b0);
 
383
        set2->vol = (set2->r1-set2->r0) * (set2->g1-set2->g0) * (set2->b1-set2->b0);
 
384
 
 
385
        return 1;
 
386
}
 
387
 
 
388
 
 
389
void Mark(struct Box *cube, int label, unsigned char *tag)
 
390
{
 
391
        ILint r, g, b;
 
392
 
 
393
        for (r = cube->r0 + 1; r <= cube->r1; r++) {
 
394
                for (g = cube->g0 + 1; g <= cube->g1; g++) {
 
395
                        for (b = cube->b0 + 1; b <= cube->b1; b++) {
 
396
                                tag[(r<<10) + (r<<6) + r + (g<<5) + g + b] = label;
 
397
                        }
 
398
                }
 
399
        }
 
400
        return;
 
401
}
 
402
 
 
403
 
 
404
ILimage *iQuantizeImage(ILimage *Image, ILuint NumCols)
 
405
{
 
406
        Box             cube[MAXCOLOR];
 
407
        ILubyte *tag = NULL;
 
408
        ILubyte lut_r[MAXCOLOR], lut_g[MAXCOLOR], lut_b[MAXCOLOR];
 
409
        ILint   next;
 
410
        ILint   weight;
 
411
        ILuint  k;
 
412
        ILfloat vv[MAXCOLOR], temp;
 
413
        //ILint color_num;
 
414
        ILubyte *NewData = NULL, *Palette = NULL;
 
415
        ILimage *TempImage = NULL, *NewImage = NULL;
 
416
        ILubyte *Ir = NULL, *Ig = NULL, *Ib = NULL;
 
417
 
 
418
        ILint num_alloced_colors; // number of colors we allocated space for in palette, as NumCols but eill not be less than 256
 
419
 
 
420
        num_alloced_colors=NumCols;
 
421
        if(num_alloced_colors<256) { num_alloced_colors=256; }
 
422
 
 
423
 
 
424
        NewImage = iCurImage;
 
425
        iCurImage = Image;
 
426
        TempImage = iConvertImage(iCurImage, IL_RGB, IL_UNSIGNED_BYTE);
 
427
        iCurImage = NewImage;
 
428
 
 
429
 
 
430
 
 
431
        if (TempImage == NULL)
 
432
                return NULL;
 
433
 
 
434
        buffer = Image->Data;
 
435
        WindW = Width = Image->Width;
 
436
        WindH = Height = Image->Height;
 
437
        WindD = Depth = Image->Depth;
 
438
        Comp = Image->Bpp;
 
439
        Qadd = NULL;
 
440
 
 
441
        //color_num = ImagePrecalculate(Image);
 
442
 
 
443
        NewData = (ILubyte*)ialloc(Image->Width * Image->Height * Image->Depth);
 
444
        Palette = (ILubyte*)ialloc(3 * num_alloced_colors);
 
445
        if (!NewData || !Palette) {
 
446
                ifree(NewData);
 
447
                ifree(Palette);
 
448
                return NULL;
 
449
        }
 
450
 
 
451
        Ir = ialloc(Width * Height * Depth);
 
452
        Ig = ialloc(Width * Height * Depth);
 
453
        Ib = ialloc(Width * Height * Depth);
 
454
        if (!Ir || !Ig || !Ib) {
 
455
                ifree(Ir);
 
456
                ifree(Ig);
 
457
                ifree(Ib);
 
458
                ifree(NewData);
 
459
                ifree(Palette);
 
460
                return NULL;
 
461
        }
 
462
 
 
463
        size = Width * Height * Depth;
 
464
        
 
465
        #ifdef ALTIVEC
 
466
            register ILuint v_size = size>>4;
 
467
            register ILuint pos = 0;
 
468
            v_size = v_size /3;
 
469
            register vector unsigned char d0,d1,d2;
 
470
            register vector unsigned char red[3],blu[3],green[3];
 
471
            
 
472
            register union{
 
473
                vector unsigned char vec;
 
474
                vector unsigned int load;
 
475
            } mask_1, mask_2, mask_3;
 
476
            
 
477
            mask_1.load = (vector unsigned int){0xFF0000FF,0x0000FF00,0x00FF0000,0xFF0000FF};
 
478
            mask_2.load = (vector unsigned int){0x00FF0000,0xFF0000FF,0x0000FF00,0x00FF0000};
 
479
            mask_2.load = (vector unsigned int){0x0000FF00,0x00FF0000,0xFF0000FF,0x0000FF00};
 
480
            
 
481
            while( v_size >= 0 ) {
 
482
                d0 = vec_ld(pos,TempImage->Data);
 
483
                d1 = vec_ld(pos+16,TempImage->Data);
 
484
                d2 = vec_ld(pos+32,TempImage->Data);
 
485
                
 
486
                red[0] =   vec_and(d0,mask_1.vec);
 
487
                green[0] = vec_and(d0,mask_2.vec);
 
488
                blu[0] =   vec_and(d0,mask_3.vec);
 
489
                
 
490
                red[1] =   vec_and(d1,mask_3.vec);
 
491
                green[1] = vec_and(d1,mask_1.vec);
 
492
                blu[1] =   vec_and(d1,mask_2.vec);
 
493
                
 
494
                red[2] =   vec_and(d2,mask_2.vec);
 
495
                green[2] = vec_and(d2,mask_3.vec);
 
496
                blu[2] =   vec_and(d2,mask_1.vec);
 
497
                
 
498
                vec_st(red[0],pos,Ir);
 
499
                vec_st(red[1],pos+16,Ir);
 
500
                vec_st(red[2],pos+32,Ir);
 
501
                
 
502
                vec_st(blu[0],pos,Ib);
 
503
                vec_st(blu[1],pos+16,Ib);
 
504
                vec_st(blu[2],pos+32,Ib);
 
505
                
 
506
                vec_st(green[0],pos,Ig);
 
507
                vec_st(green[1],pos+16,Ig);
 
508
                vec_st(green[2],pos+32,Ig);
 
509
                
 
510
                pos += 48;
 
511
            }
 
512
            size -= pos;
 
513
        #endif
 
514
 
 
515
        for (k = 0; k < size; k++) {
 
516
                Ir[k] = TempImage->Data[k * 3];
 
517
                Ig[k] = TempImage->Data[k * 3 + 1];
 
518
                Ib[k] = TempImage->Data[k * 3 + 2];
 
519
        }
 
520
        
 
521
        #ifdef ALTIVEC
 
522
           size = Width * Height * Depth;
 
523
        #endif
 
524
        
 
525
        // Set new colors number
 
526
        K = NumCols;
 
527
 
 
528
        if (K <= 256) {
 
529
                // Begin Wu's color quantization algorithm
 
530
 
 
531
                // May have "leftovers" from a previous run.
 
532
                
 
533
                imemclear(gm2, 33 * 33 * 33 * sizeof(ILfloat));
 
534
                imemclear(wt, 33 * 33 * 33 * sizeof(ILint));
 
535
                imemclear(mr, 33 * 33 * 33 * sizeof(ILint));
 
536
                imemclear(mg, 33 * 33 * 33 * sizeof(ILint));
 
537
                imemclear(mb, 33 * 33 * 33 * sizeof(ILint));
 
538
                
 
539
                if (!Hist3d(Ir, Ig, Ib, (ILint*)wt, (ILint*)mr, (ILint*)mg, (ILint*)mb, (ILfloat*)gm2))
 
540
                        goto error_label;
 
541
 
 
542
                M3d((ILint*)wt, (ILint*)mr, (ILint*)mg, (ILint*)mb, (ILfloat*)gm2);
 
543
 
 
544
                cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
 
545
                cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
 
546
                next = 0;
 
547
                for (i = 1; i < K; ++i) {
 
548
                        if (Cut(&cube[next], &cube[i])) { // volume test ensures we won't try to cut one-cell Box */
 
549
                                vv[next] = (cube[next].vol>1) ? Var(&cube[next]) : 0.0f;
 
550
                                vv[i] = (cube[i].vol>1) ? Var(&cube[i]) : 0.0f;
 
551
                        }
 
552
                        else {
 
553
                                vv[next] = 0.0;   // don't try to split this Box again
 
554
                                i--;              // didn't create Box i
 
555
                        }
 
556
                        next = 0;
 
557
                        temp = vv[0];
 
558
                        for (k = 1; (ILint)k <= i; ++k) {
 
559
                                if (vv[k] > temp) {
 
560
                                        temp = vv[k]; next = k;
 
561
                                }
 
562
                        }
 
563
                                
 
564
                        if (temp <= 0.0) {
 
565
                                K = i+1;
 
566
                                // Only got K Boxes
 
567
                                break;
 
568
                        }
 
569
                }
 
570
 
 
571
                tag = (ILubyte*)ialloc(33 * 33 * 33 * sizeof(ILubyte));
 
572
                if (tag == NULL)
 
573
                        goto error_label;
 
574
                for (k = 0; (ILint)k < K; k++) {
 
575
                        Mark(&cube[k], k, tag);
 
576
                        weight = Vol(&cube[k], wt);
 
577
                        if (weight) {
 
578
                                lut_r[k] = (ILubyte)(Vol(&cube[k], mr) / weight);
 
579
                                lut_g[k] = (ILubyte)(Vol(&cube[k], mg) / weight);
 
580
                                lut_b[k] = (ILubyte)(Vol(&cube[k], mb) / weight);
 
581
                        }
 
582
                        else {
 
583
                                // Bogus Box
 
584
                                lut_r[k] = lut_g[k] = lut_b[k] = 0;             
 
585
                        }
 
586
                }
 
587
 
 
588
                for (i = 0; i < (ILint)size; i++) {
 
589
                        NewData[i] = tag[Qadd[i]];
 
590
                }
 
591
                ifree(tag);
 
592
                ifree(Qadd);
 
593
 
 
594
                for (k = 0; k < NumCols; k++) {
 
595
                        Palette[k * 3]     = lut_b[k];
 
596
                        Palette[k * 3 + 1] = lut_g[k];
 
597
                        Palette[k * 3 + 2] = lut_r[k];
 
598
                }
 
599
        }
 
600
        else { // If colors more than 256
 
601
                // Begin Octree quantization
 
602
                //Quant_Octree(Image->Width, Image->Height, Image->Bpp, buffer2, NewData, Palette, K);
 
603
                ilSetError(IL_INTERNAL_ERROR);  // Can't get much more specific than this.
 
604
                goto error_label;
 
605
        }
 
606
 
 
607
        ifree(Ig);
 
608
        ifree(Ib);
 
609
        ifree(Ir);
 
610
        ilCloseImage(TempImage);
 
611
 
 
612
        NewImage = (ILimage*)icalloc(sizeof(ILimage), 1);
 
613
        if (NewImage == NULL) {
 
614
                return NULL;
 
615
        }
 
616
        ilCopyImageAttr(NewImage, Image);
 
617
        NewImage->Bpp = 1;
 
618
        NewImage->Bps = Image->Width;
 
619
        NewImage->SizeOfPlane = NewImage->Bps * Image->Height;
 
620
        NewImage->SizeOfData = NewImage->SizeOfPlane;
 
621
        NewImage->Format = IL_COLOUR_INDEX;
 
622
        NewImage->Type = IL_UNSIGNED_BYTE;
 
623
 
 
624
        NewImage->Pal.Palette = Palette;
 
625
        NewImage->Pal.PalSize = 256 * 3;
 
626
        NewImage->Pal.PalType = IL_PAL_BGR24;
 
627
        NewImage->Data = NewData;
 
628
 
 
629
        return NewImage;
 
630
 
 
631
error_label:
 
632
        ifree(NewData);
 
633
        ifree(Palette);
 
634
        ifree(Ig);
 
635
        ifree(Ib);
 
636
        ifree(Ir);
 
637
        ifree(tag);
 
638
        ifree(Qadd);
 
639
        return NULL;
 
640
}