1 /* median.c: median cut - reducing a high color bitmap to certain number of colors
2
3 Copyright (C) 2001, 2002 Martin Weber
4
5 This library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public License
7 as published by the Free Software Foundation; either version 2.1 of
8 the License, or (at your option) any later version.
9
10 This library is distributed in the hope that it will be useful, but
11 WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with this library; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
18 USA. */
19
20 #ifdef HAVE_CONFIG_H
21 #include "config.h"
22 #endif /* Def: HAVE_CONFIG_H */
23
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include "logreport.h"
28 #include "xstd.h"
29 #include "quantize.h"
30
31 #define MAXNUMCOLORS 256
32
33 #if 0
34 #define R_SCALE
35 #define G_SCALE
36 #define B_SCALE
37 #else
38
39 /* scale RGB distances by *2,*3,*1 */
40 #define R_SCALE <<1
41 #define G_SCALE *3
42 #define B_SCALE
43 #endif
44
45 #define BITS_IN_SAMPLE 8
46
47 #define R_SHIFT (BITS_IN_SAMPLE - PRECISION_R)
48 #define G_SHIFT (BITS_IN_SAMPLE - PRECISION_G)
49 #define B_SHIFT (BITS_IN_SAMPLE - PRECISION_B)
50
51 typedef struct {
52 /* The bounds of the box (inclusive); expressed as histogram indexes */
53 int Rmin, Rmax;
54 int Gmin, Gmax;
55 int Bmin, Bmax;
56 /* The volume (actually 2-norm) of the box */
57 int volume;
58 /* The number of nonzero histogram cells within this box */
59 long colorcount;
60 } box, *boxptr;
61
zero_histogram_rgb(Histogram histogram)62 static void zero_histogram_rgb(Histogram histogram)
63 {
64 int r, g, b;
65 for (r = 0; r < HIST_R_ELEMS; r++)
66 for (g = 0; g < HIST_G_ELEMS; g++)
67 for (b = 0; b < HIST_B_ELEMS; b++)
68 histogram[r * MR + g * MG + b] = 0;
69 }
70
generate_histogram_rgb(Histogram histogram,at_bitmap * image,const at_color * ignoreColor)71 static void generate_histogram_rgb(Histogram histogram, at_bitmap * image, const at_color * ignoreColor)
72 {
73 unsigned char *src = image->bitmap;
74 int num_elems;
75 ColorFreq *col;
76
77 num_elems = AT_BITMAP_WIDTH(image) * AT_BITMAP_HEIGHT(image);
78 zero_histogram_rgb(histogram);
79
80 switch (AT_BITMAP_PLANES(image)) {
81 case 3:
82 while (num_elems--) {
83 /* If we have an ignorecolor, skip it. */
84 if (ignoreColor) {
85 if ((src[0] == ignoreColor->r)
86 && (src[1] == ignoreColor->g)
87 && (src[2] == ignoreColor->b)) {
88 src += 3;
89 continue;
90 }
91 }
92 col = &histogram[(src[0] >> R_SHIFT) * MR + (src[1] >> G_SHIFT) * MG + (src[2] >> B_SHIFT)];
93 (*col)++;
94 src += 3;
95 }
96 break;
97
98 case 1:
99 while (--num_elems >= 0) {
100 if (ignoreColor && src[num_elems] == ignoreColor->r)
101 continue;
102 col = &histogram[(src[num_elems] >> R_SHIFT) * MR + (src[num_elems] >> G_SHIFT) * MG + (src[num_elems] >> B_SHIFT)];
103 (*col)++;
104 }
105 break;
106 default:
107 /* To avoid compiler warning */ ;
108 }
109 }
110
find_biggest_volume(boxptr boxlist,int numboxes)111 static boxptr find_biggest_volume(boxptr boxlist, int numboxes)
112 /* Find the splittable box with the largest (scaled) volume */
113 /* Returns 0 if no splittable boxes remain */
114 {
115 boxptr boxp;
116 int i;
117 int maxv = 0;
118 boxptr which = 0;
119
120 for (i = 0, boxp = boxlist; i < numboxes; i++, boxp++) {
121 if (boxp->volume > maxv) {
122 which = boxp;
123 maxv = boxp->volume;
124 }
125 }
126
127 return which;
128 }
129
update_box_rgb(Histogram histogram,boxptr boxp)130 static void update_box_rgb(Histogram histogram, boxptr boxp)
131 /* Shrink the min/max bounds of a box to enclose only nonzero elements, */
132 /* and recompute its volume and population */
133 {
134 ColorFreq *histp;
135 int R, G, B;
136 int Rmin, Rmax, Gmin, Gmax, Bmin, Bmax;
137 int dist0, dist1, dist2;
138 long ccount;
139
140 Rmin = boxp->Rmin;
141 Rmax = boxp->Rmax;
142 Gmin = boxp->Gmin;
143 Gmax = boxp->Gmax;
144 Bmin = boxp->Bmin;
145 Bmax = boxp->Bmax;
146
147 if (Rmax > Rmin)
148 for (R = Rmin; R <= Rmax; R++)
149 for (G = Gmin; G <= Gmax; G++) {
150 histp = histogram + R * MR + G * MG + Bmin;
151 for (B = Bmin; B <= Bmax; B++)
152 if (*histp++ != 0) {
153 boxp->Rmin = Rmin = R;
154 goto have_Rmin;
155 }
156 }
157 have_Rmin:
158 if (Rmax > Rmin)
159 for (R = Rmax; R >= Rmin; R--)
160 for (G = Gmin; G <= Gmax; G++) {
161 histp = histogram + R * MR + G * MG + Bmin;
162 for (B = Bmin; B <= Bmax; B++)
163 if (*histp++ != 0) {
164 boxp->Rmax = Rmax = R;
165 goto have_Rmax;
166 }
167 }
168 have_Rmax:
169 if (Gmax > Gmin)
170 for (G = Gmin; G <= Gmax; G++)
171 for (R = Rmin; R <= Rmax; R++) {
172 histp = histogram + R * MR + G * MG + Bmin;
173 for (B = Bmin; B <= Bmax; B++)
174 if (*histp++ != 0) {
175 boxp->Gmin = Gmin = G;
176 goto have_Gmin;
177 }
178 }
179 have_Gmin:
180 if (Gmax > Gmin)
181 for (G = Gmax; G >= Gmin; G--)
182 for (R = Rmin; R <= Rmax; R++) {
183 histp = histogram + R * MR + G * MG + Bmin;
184 for (B = Bmin; B <= Bmax; B++)
185 if (*histp++ != 0) {
186 boxp->Gmax = Gmax = G;
187 goto have_Gmax;
188 }
189 }
190 have_Gmax:
191 if (Bmax > Bmin)
192 for (B = Bmin; B <= Bmax; B++)
193 for (R = Rmin; R <= Rmax; R++) {
194 histp = histogram + R * MR + Gmin * MG + B;
195 for (G = Gmin; G <= Gmax; G++, histp += MG)
196 if (*histp != 0) {
197 boxp->Bmin = Bmin = B;
198 goto have_Bmin;
199 }
200 }
201 have_Bmin:
202 if (Bmax > Bmin)
203 for (B = Bmax; B >= Bmin; B--)
204 for (R = Rmin; R <= Rmax; R++) {
205 histp = histogram + R * MR + Gmin * MG + B;
206 for (G = Gmin; G <= Gmax; G++, histp += MG)
207 if (*histp != 0) {
208 boxp->Bmax = Bmax = B;
209 goto have_Bmax;
210 }
211 }
212 have_Bmax:
213
214 /* Update box volume.
215 * We use 2-norm rather than real volume here; this biases the method
216 * against making long narrow boxes, and it has the side benefit that
217 * a box is splittable iff norm > 0.
218 * Since the differences are expressed in histogram-cell units,
219 * we have to shift back to JSAMPLE units to get consistent distances;
220 * after which, we scale according to the selected distance scale factors.
221 */
222 dist0 = Rmax - Rmin;
223 dist1 = Gmax - Gmin;
224 dist2 = Bmax - Bmin;
225 boxp->volume = dist0 * dist0 + dist1 * dist1 + dist2 * dist2;
226
227 /* Now scan remaining volume of box and compute population */
228 ccount = 0;
229 for (R = Rmin; R <= Rmax; R++)
230 for (G = Gmin; G <= Gmax; G++) {
231 histp = histogram + R * MR + G * MG + Bmin;
232 for (B = Bmin; B <= Bmax; B++, histp++)
233 if (*histp != 0) {
234 ccount++;
235 }
236 }
237
238 boxp->colorcount = ccount;
239 }
240
median_cut_rgb(Histogram histogram,boxptr boxlist,int numboxes,int desired_colors)241 static int median_cut_rgb(Histogram histogram, boxptr boxlist, int numboxes, int desired_colors)
242 /* Repeatedly select and split the largest box until we have enough boxes */
243 {
244 int n, lb;
245 int R, G, B, cmax;
246 boxptr b1, b2;
247
248 while (numboxes < desired_colors) {
249 /* Select box to split.
250 * Current algorithm: by population for first half, then by volume.
251 */
252 b1 = find_biggest_volume(boxlist, numboxes);
253
254 if (b1 == 0) /* no splittable boxes left! */
255 break;
256 b2 = boxlist + numboxes; /* where new box will go */
257 /* Copy the color bounds to the new box. */
258 b2->Rmax = b1->Rmax;
259 b2->Gmax = b1->Gmax;
260 b2->Bmax = b1->Bmax;
261 b2->Rmin = b1->Rmin;
262 b2->Gmin = b1->Gmin;
263 b2->Bmin = b1->Bmin;
264 /* Choose which axis to split the box on.
265 * Current algorithm: longest scaled axis.
266 * See notes in update_box about scaling distances.
267 */
268 R = b1->Rmax - b1->Rmin;
269 G = b1->Gmax - b1->Gmin;
270 B = b1->Bmax - b1->Bmin;
271 /* We want to break any ties in favor of green, then red, blue last.
272 */
273 cmax = G;
274 n = 1;
275 if (R > cmax) {
276 cmax = R;
277 n = 0;
278 }
279 if (B > cmax) {
280 n = 2;
281 }
282 /* Choose split point along selected axis, and update box bounds.
283 * Current algorithm: split at halfway point.
284 * (Since the box has been shrunk to minimum volume,
285 * any split will produce two nonempty subboxes.)
286 * Note that lb value is max for lower box, so must be < old max.
287 */
288 switch (n) {
289 case 0:
290 lb = (b1->Rmax + b1->Rmin) / 2;
291 b1->Rmax = lb;
292 b2->Rmin = lb + 1;
293 break;
294 case 1:
295 lb = (b1->Gmax + b1->Gmin) / 2;
296 b1->Gmax = lb;
297 b2->Gmin = lb + 1;
298 break;
299 case 2:
300 lb = (b1->Bmax + b1->Bmin) / 2;
301 b1->Bmax = lb;
302 b2->Bmin = lb + 1;
303 break;
304 }
305 /* Update stats for boxes */
306 update_box_rgb(histogram, b1);
307 update_box_rgb(histogram, b2);
308 numboxes++;
309 }
310 return numboxes;
311 }
312
compute_color_rgb(QuantizeObj * quantobj,Histogram histogram,boxptr boxp,int icolor)313 static void compute_color_rgb(QuantizeObj * quantobj, Histogram histogram, boxptr boxp, int icolor)
314 /* Compute representative color for a box, put it in colormap[icolor] */
315 {
316 /* Current algorithm: mean weighted by pixels (not colors) */
317 /* Note it is important to get the rounding correct! */
318 ColorFreq *histp;
319 int R, G, B;
320 int Rmin, Rmax;
321 int Gmin, Gmax;
322 int Bmin, Bmax;
323 unsigned long count;
324 unsigned long total = 0;
325 unsigned long Rtotal = 0;
326 unsigned long Gtotal = 0;
327 unsigned long Btotal = 0;
328
329 Rmin = boxp->Rmin;
330 Rmax = boxp->Rmax;
331 Gmin = boxp->Gmin;
332 Gmax = boxp->Gmax;
333 Bmin = boxp->Bmin;
334 Bmax = boxp->Bmax;
335
336 for (R = Rmin; R <= Rmax; R++)
337 for (G = Gmin; G <= Gmax; G++) {
338 histp = histogram + R * MR + G * MG + Bmin;
339 for (B = Bmin; B <= Bmax; B++) {
340 if ((count = *histp++) != 0) {
341 total += count;
342 Rtotal += ((R << R_SHIFT) + ((1 << R_SHIFT) >> 1)) * count;
343 Gtotal += ((G << G_SHIFT) + ((1 << G_SHIFT) >> 1)) * count;
344 Btotal += ((B << B_SHIFT) + ((1 << B_SHIFT) >> 1)) * count;
345 }
346 }
347 }
348
349 quantobj->cmap[icolor].r = (unsigned char)((Rtotal + (total >> 1)) / total);
350 quantobj->cmap[icolor].g = (unsigned char)((Gtotal + (total >> 1)) / total);
351 quantobj->cmap[icolor].b = (unsigned char)((Btotal + (total >> 1)) / total);
352 quantobj->freq[icolor] = total;
353 }
354
select_colors_rgb(QuantizeObj * quantobj,Histogram histogram)355 static void select_colors_rgb(QuantizeObj * quantobj, Histogram histogram)
356 /* Master routine for color selection */
357 {
358 boxptr boxlist;
359 int numboxes;
360 int desired = quantobj->desired_number_of_colors;
361 int i;
362
363 /* Allocate workspace for box list */
364 XMALLOC(boxlist, desired * sizeof(box));
365
366 /* Initialize one box containing whole space */
367 numboxes = 1;
368 boxlist[0].Rmin = 0;
369 boxlist[0].Rmax = (1 << PRECISION_R) - 1;
370 boxlist[0].Gmin = 0;
371 boxlist[0].Gmax = (1 << PRECISION_G) - 1;
372 boxlist[0].Bmin = 0;
373 boxlist[0].Bmax = (1 << PRECISION_B) - 1;
374 /* Shrink it to actually-used volume and set its statistics */
375 update_box_rgb(histogram, boxlist);
376 /* Perform median-cut to produce final box list */
377 numboxes = median_cut_rgb(histogram, boxlist, numboxes, desired);
378 quantobj->actual_number_of_colors = numboxes;
379 /* Compute the representative color for each box, fill colormap */
380 for (i = 0; i < numboxes; i++)
381 compute_color_rgb(quantobj, histogram, boxlist + i, i);
382 free(boxlist);
383 }
384
385 /*
386 * These routines are concerned with the time-critical task of mapping input
387 * colors to the nearest color in the selected colormap.
388 *
389 * We re-use the histogram space as an "inverse color map", essentially a
390 * cache for the results of nearest-color searches. All colors within a
391 * histogram cell will be mapped to the same colormap entry, namely the one
392 * closest to the cell's center. This may not be quite the closest entry to
393 * the actual input color, but it's almost as good. A zero in the cache
394 * indicates we haven't found the nearest color for that cell yet; the array
395 * is cleared to zeroes before starting the mapping pass. When we find the
396 * nearest color for a cell, its colormap index plus one is recorded in the
397 * cache for future use. The pass2 scanning routines call fill_inverse_cmap
398 * when they need to use an unfilled entry in the cache.
399 *
400 * Our method of efficiently finding nearest colors is based on the "locally
401 * sorted search" idea described by Heckbert and on the incremental distance
402 * calculation described by Spencer W. Thomas in chapter III.1 of Graphics
403 * Gems II (James Arvo, ed. Academic Press, 1991). Thomas points out that
404 * the distances from a given colormap entry to each cell of the histogram can
405 * be computed quickly using an incremental method: the differences between
406 * distances to adjacent cells themselves differ by a constant. This allows a
407 * fairly fast implementation of the "brute force" approach of computing the
408 * distance from every colormap entry to every histogram cell. Unfortunately,
409 * it needs a work array to hold the best-distance-so-far for each histogram
410 * cell (because the inner loop has to be over cells, not colormap entries).
411 * The work array elements have to be ints, so the work array would need
412 * 256Kb at our recommended precision. This is not feasible in DOS machines.
413
414 [ 256*1024/4 = 65,536 ]
415
416 * To get around these problems, we apply Thomas' method to compute the
417 * nearest colors for only the cells within a small subbox of the histogram.
418 * The work array need be only as big as the subbox, so the memory usage
419 * problem is solved. Furthermore, we need not fill subboxes that are never
420 * referenced in pass2; many images use only part of the color gamut, so a
421 * fair amount of work is saved. An additional advantage of this
422 * approach is that we can apply Heckbert's locality criterion to quickly
423 * eliminate colormap entries that are far away from the subbox; typically
424 * three-fourths of the colormap entries are rejected by Heckbert's criterion,
425 * and we need not compute their distances to individual cells in the subbox.
426 * The speed of this approach is heavily influenced by the subbox size: too
427 * small means too much overhead, too big loses because Heckbert's criterion
428 * can't eliminate as many colormap entries. Empirically the best subbox
429 * size seems to be about 1/512th of the histogram (1/8th in each direction).
430 *
431 * Thomas' article also describes a refined method which is asymptotically
432 * faster than the brute-force method, but it is also far more complex and
433 * cannot efficiently be applied to small subboxes. It is therefore not
434 * useful for programs intended to be portable to DOS machines. On machines
435 * with plenty of memory, filling the whole histogram in one shot with Thomas'
436 * refined method might be faster than the present code --- but then again,
437 * it might not be any faster, and it's certainly more complicated.
438 */
439
440 /* log2(histogram cells in update box) for each axis; this can be adjusted */
441 #define BOX_R_LOG (PRECISION_R-3)
442 #define BOX_G_LOG (PRECISION_G-3)
443 #define BOX_B_LOG (PRECISION_B-3)
444
445 #define BOX_R_ELEMS (1<<BOX_R_LOG) /* # of hist cells in update box */
446 #define BOX_G_ELEMS (1<<BOX_G_LOG)
447 #define BOX_B_ELEMS (1<<BOX_B_LOG)
448
449 #define BOX_R_SHIFT (R_SHIFT + BOX_R_LOG)
450 #define BOX_G_SHIFT (G_SHIFT + BOX_G_LOG)
451 #define BOX_B_SHIFT (B_SHIFT + BOX_B_LOG)
452
453 /*
454 * The next three routines implement inverse colormap filling. They could
455 * all be folded into one big routine, but splitting them up this way saves
456 * some stack space (the mindist[] and bestdist[] arrays need not coexist)
457 * and may allow some compilers to produce better code by registerizing more
458 * inner-loop variables.
459 */
460
find_nearby_colors(QuantizeObj * quantobj,int minR,int minG,int minB,int * colorlist)461 static int find_nearby_colors(QuantizeObj * quantobj, int minR, int minG, int minB, int *colorlist)
462 /* Locate the colormap entries close enough to an update box to be candidates
463 * for the nearest entry to some cell(s) in the update box. The update box
464 * is specified by the center coordinates of its first cell. The number of
465 * candidate colormap entries is returned, and their colormap indexes are
466 * placed in colorlist[].
467 * This routine uses Heckbert's "locally sorted search" criterion to select
468 * the colors that need further consideration.
469 */
470 {
471 int numcolors = quantobj->actual_number_of_colors;
472 int maxR, maxG, maxB;
473 int centerR, centerG, centerB;
474 int i, x, ncolors;
475 int minmaxdist, min_dist = 0, max_dist, tdist;
476 int mindist[MAXNUMCOLORS]; /* min distance to colormap entry i */
477
478 /* Compute TRUE coordinates of update box's upper corner and center.
479 * Actually we compute the coordinates of the center of the upper-corner
480 * histogram cell, which are the upper bounds of the volume we care about.
481 * Note that since ">>" rounds down, the "center" values may be closer to
482 * min than to max; hence comparisons to them must be "<=", not "<".
483 */
484 maxR = minR + ((1 << BOX_R_SHIFT) - (1 << R_SHIFT));
485 centerR = (minR + maxR) >> 1;
486 maxG = minG + ((1 << BOX_G_SHIFT) - (1 << G_SHIFT));
487 centerG = (minG + maxG) >> 1;
488 maxB = minB + ((1 << BOX_B_SHIFT) - (1 << B_SHIFT));
489 centerB = (minB + maxB) >> 1;
490
491 /* For each color in colormap, find:
492 * 1. its minimum squared-distance to any point in the update box
493 * (zero if color is within update box);
494 * 2. its maximum squared-distance to any point in the update box.
495 * Both of these can be found by considering only the corners of the box.
496 * We save the minimum distance for each color in mindist[];
497 * only the smallest maximum distance is of interest.
498 */
499 minmaxdist = 0x7FFFFFFFL;
500
501 for (i = 0; i < numcolors; i++) {
502 /* We compute the squared-R-distance term, then add in the other two. */
503 x = quantobj->cmap[i].r;
504 if (x < minR) {
505 tdist = (x - minR) R_SCALE;
506 min_dist = tdist * tdist;
507 tdist = (x - maxR) R_SCALE;
508 max_dist = tdist * tdist;
509 } else if (x > maxR) {
510 tdist = (x - maxR) R_SCALE;
511 min_dist = tdist * tdist;
512 tdist = (x - minR) R_SCALE;
513 max_dist = tdist * tdist;
514 } else {
515 /* within cell range so no contribution to min_dist */
516 min_dist = 0;
517 if (x <= centerR) {
518 tdist = (x - maxR) R_SCALE;
519 max_dist = tdist * tdist;
520 } else {
521 tdist = (x - minR) R_SCALE;
522 max_dist = tdist * tdist;
523 }
524 }
525
526 x = quantobj->cmap[i].g;
527 if (x < minG) {
528 tdist = (x - minG) G_SCALE;
529 min_dist += tdist * tdist;
530 tdist = (x - maxG) G_SCALE;
531 max_dist += tdist * tdist;
532 } else if (x > maxG) {
533 tdist = (x - maxG) G_SCALE;
534 min_dist += tdist * tdist;
535 tdist = (x - minG) G_SCALE;
536 max_dist += tdist * tdist;
537 } else {
538 /* within cell range so no contribution to min_dist */
539 if (x <= centerG) {
540 tdist = (x - maxG) G_SCALE;
541 max_dist += tdist * tdist;
542 } else {
543 tdist = (x - minG) G_SCALE;
544 max_dist += tdist * tdist;
545 }
546 }
547
548 x = quantobj->cmap[i].b;
549 if (x < minB) {
550 tdist = (x - minB) B_SCALE;
551 min_dist += tdist * tdist;
552 tdist = (x - maxB) B_SCALE;
553 max_dist += tdist * tdist;
554 } else if (x > maxB) {
555 tdist = (x - maxB) B_SCALE;
556 min_dist += tdist * tdist;
557 tdist = (x - minB) B_SCALE;
558 max_dist += tdist * tdist;
559 } else {
560 /* within cell range so no contribution to min_dist */
561 if (x <= centerB) {
562 tdist = (x - maxB) B_SCALE;
563 max_dist += tdist * tdist;
564 } else {
565 tdist = (x - minB) B_SCALE;
566 max_dist += tdist * tdist;
567 }
568 }
569
570 mindist[i] = min_dist; /* save away the results */
571 if (max_dist < minmaxdist)
572 minmaxdist = max_dist;
573 }
574
575 /* Now we know that no cell in the update box is more than minmaxdist
576 * away from some colormap entry. Therefore, only colors that are
577 * within minmaxdist of some part of the box need be considered.
578 */
579 ncolors = 0;
580 for (i = 0; i < numcolors; i++) {
581 if (mindist[i] <= minmaxdist)
582 colorlist[ncolors++] = i;
583 }
584 return ncolors;
585 }
586
find_best_colors(QuantizeObj * quantobj,int minR,int minG,int minB,int numcolors,int * colorlist,int * bestcolor)587 static void find_best_colors(QuantizeObj * quantobj, int minR, int minG, int minB, int numcolors, int *colorlist, int *bestcolor)
588 /* Find the closest colormap entry for each cell in the update box,
589 given the list of candidate colors prepared by find_nearby_colors.
590 Return the indexes of the closest entries in the bestcolor[] array.
591 This routine uses Thomas' incremental distance calculation method to
592 find the distance from a colormap entry to successive cells in the box.
593 */
594 {
595 int iR, iG, iB;
596 int i, icolor;
597 int *bptr; /* pointer into bestdist[] array */
598 int *cptr; /* pointer into bestcolor[] array */
599 int dist0, dist1; /* initial distance values */
600 int dist2; /* current distance in inner loop */
601 int xx0, xx1; /* distance increments */
602 int xx2;
603 int inR, inG, inB; /* initial values for increments */
604
605 /* This array holds the distance to the nearest-so-far color for each cell */
606 int bestdist[BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS];
607
608 /* Initialize best-distance for each cell of the update box */
609 bptr = bestdist;
610 for (i = BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS - 1; i >= 0; i--)
611 *bptr++ = 0x7FFFFFFFL;
612
613 /* For each color selected by find_nearby_colors,
614 * compute its distance to the center of each cell in the box.
615 * If that's less than best-so-far, update best distance and color number.
616 */
617
618 /* Nominal steps between cell centers ("x" in Thomas article) */
619 #define STEP_R ((1 << R_SHIFT) R_SCALE)
620 #define STEP_G ((1 << G_SHIFT) G_SCALE)
621 #define STEP_B ((1 << B_SHIFT) B_SCALE)
622
623 for (i = 0; i < numcolors; i++) {
624 icolor = colorlist[i];
625 /* Compute (square of) distance from minR/G/B to this color */
626 inR = (minR - quantobj->cmap[icolor].r) R_SCALE;
627 dist0 = inR * inR;
628 inG = (minG - quantobj->cmap[icolor].g) G_SCALE;
629 dist0 += inG * inG;
630 inB = (minB - quantobj->cmap[icolor].b) B_SCALE;
631 dist0 += inB * inB;
632 /* Form the initial difference increments */
633 inR = inR * (2 * STEP_R) + STEP_R * STEP_R;
634 inG = inG * (2 * STEP_G) + STEP_G * STEP_G;
635 inB = inB * (2 * STEP_B) + STEP_B * STEP_B;
636 /* Now loop over all cells in box, updating distance per Thomas method */
637 bptr = bestdist;
638 cptr = bestcolor;
639 xx0 = inR;
640 for (iR = BOX_R_ELEMS - 1; iR >= 0; iR--) {
641 dist1 = dist0;
642 xx1 = inG;
643 for (iG = BOX_G_ELEMS - 1; iG >= 0; iG--) {
644 dist2 = dist1;
645 xx2 = inB;
646 for (iB = BOX_B_ELEMS - 1; iB >= 0; iB--) {
647 if (dist2 < *bptr) {
648 *bptr = dist2;
649 *cptr = icolor;
650 }
651 dist2 += xx2;
652 xx2 += 2 * STEP_B * STEP_B;
653 bptr++;
654 cptr++;
655 }
656 dist1 += xx1;
657 xx1 += 2 * STEP_G * STEP_G;
658 }
659 dist0 += xx0;
660 xx0 += 2 * STEP_R * STEP_R;
661 }
662 }
663 }
664
fill_inverse_cmap_rgb(QuantizeObj * quantobj,Histogram histogram,int R,int G,int B)665 static void fill_inverse_cmap_rgb(QuantizeObj * quantobj, Histogram histogram, int R, int G, int B)
666 /* Fill the inverse-colormap entries in the update box that contains
667 histogram cell R/G/B. (Only that one cell MUST be filled, but
668 we can fill as many others as we wish.) */
669 {
670 int minR, minG, minB; /* lower left corner of update box */
671 int iR, iG, iB;
672 int *cptr; /* pointer into bestcolor[] array */
673 ColorFreq *cachep; /* pointer into main cache array */
674 /* This array lists the candidate colormap indexes. */
675 int colorlist[MAXNUMCOLORS];
676 int numcolors; /* number of candidate colors */
677 /* This array holds the actually closest colormap index for each cell. */
678 int bestcolor[BOX_R_ELEMS * BOX_G_ELEMS * BOX_B_ELEMS];
679
680 /* Convert cell coordinates to update box ID */
681 R >>= BOX_R_LOG;
682 G >>= BOX_G_LOG;
683 B >>= BOX_B_LOG;
684
685 /* Compute TRUE coordinates of update box's origin corner.
686 * Actually we compute the coordinates of the center of the corner
687 * histogram cell, which are the lower bounds of the volume we care about.
688 */
689 minR = (R << BOX_R_SHIFT) + ((1 << R_SHIFT) >> 1);
690 minG = (G << BOX_G_SHIFT) + ((1 << G_SHIFT) >> 1);
691 minB = (B << BOX_B_SHIFT) + ((1 << B_SHIFT) >> 1);
692
693 /* Determine which colormap entries are close enough to be candidates
694 * for the nearest entry to some cell in the update box.
695 */
696 numcolors = find_nearby_colors(quantobj, minR, minG, minB, colorlist);
697
698 /* Determine the actually nearest colors. */
699 find_best_colors(quantobj, minR, minG, minB, numcolors, colorlist, bestcolor);
700
701 /* Save the best color numbers (plus 1) in the main cache array */
702 R <<= BOX_R_LOG; /* convert ID back to base cell indexes */
703 G <<= BOX_G_LOG;
704 B <<= BOX_B_LOG;
705 cptr = bestcolor;
706 for (iR = 0; iR < BOX_R_ELEMS; iR++) {
707 for (iG = 0; iG < BOX_G_ELEMS; iG++) {
708 cachep = &histogram[(R + iR) * MR + (G + iG) * MG + B];
709 for (iB = 0; iB < BOX_B_ELEMS; iB++) {
710 *cachep++ = (*cptr++) + 1;
711 }
712 }
713 }
714 }
715
716 /* This is pass 1 */
median_cut_pass1_rgb(QuantizeObj * quantobj,at_bitmap * image,const at_color * ignoreColor)717 static void median_cut_pass1_rgb(QuantizeObj * quantobj, at_bitmap * image, const at_color * ignoreColor)
718 {
719 generate_histogram_rgb(quantobj->histogram, image, ignoreColor);
720 select_colors_rgb(quantobj, quantobj->histogram);
721 }
722
723 /* Map some rows of pixels to the output colormapped representation. */
median_cut_pass2_rgb(QuantizeObj * quantobj,at_bitmap * image,const at_color * bgColor)724 static void median_cut_pass2_rgb(QuantizeObj * quantobj, at_bitmap * image, const at_color * bgColor)
725 /* This version performs no dithering */
726 {
727 Histogram histogram = quantobj->histogram;
728 ColorFreq *cachep;
729 int R, G, B;
730 int origR, origG, origB;
731 int row, col;
732 int spp = AT_BITMAP_PLANES(image);
733 int width = AT_BITMAP_WIDTH(image);
734 int height = AT_BITMAP_HEIGHT(image);
735 unsigned char *src, *dest;
736 at_color bg_color = { 0xff, 0xff, 0xff };
737
738 zero_histogram_rgb(histogram);
739
740 if (bgColor) {
741 /* Find the nearest colormap entry for the background color. */
742 R = bgColor->r >> R_SHIFT;
743 G = bgColor->g >> G_SHIFT;
744 B = bgColor->b >> B_SHIFT;
745 cachep = &histogram[R * MR + G * MG + B];
746 if (*cachep == 0)
747 fill_inverse_cmap_rgb(quantobj, histogram, R, G, B);
748 bg_color = quantobj->cmap[*cachep - 1];
749 }
750
751 src = dest = image->bitmap;
752 if (spp == 3) {
753 for (row = 0; row < height; row++) {
754 for (col = 0; col < width; col++) {
755 /* get pixel value and index into the cache */
756 origR = (*src++);
757 origG = (*src++);
758 origB = (*src++);
759
760 /*
761 if (origR > 253 && origG > 253 && origB > 253)
762 {
763 (*dest++) = 255; (*dest++) = 255; (*dest++) = 255;
764 continue;
765 }
766 */
767
768 /* get pixel value and index into the cache */
769 R = origR >> R_SHIFT;
770 G = origG >> G_SHIFT;
771 B = origB >> B_SHIFT;
772 cachep = &histogram[R * MR + G * MG + B];
773 /* If we have not seen this color before, find nearest
774 colormap entry and update the cache */
775 if (*cachep == 0) {
776 fill_inverse_cmap_rgb(quantobj, histogram, R, G, B);
777 }
778 /* Now emit the colormap index for this cell */
779 dest[0] = quantobj->cmap[*cachep - 1].r;
780 dest[1] = quantobj->cmap[*cachep - 1].g;
781 dest[2] = quantobj->cmap[*cachep - 1].b;
782
783 /* If the colormap entry for this pixel is the same as the
784 background's colormap entry, set the pixel to the
785 background color. */
786 if (bgColor && (dest[0] == bg_color.r && dest[1] == bg_color.g && dest[2] == bg_color.b)) {
787 dest[0] = bgColor->r;
788 dest[1] = bgColor->g;
789 dest[2] = bgColor->b;
790 }
791 dest += 3;
792 }
793 }
794 } else if (spp == 1) {
795 long idx = width * height;
796 while (--idx >= 0) {
797 origR = src[idx];
798 R = origR >> R_SHIFT;
799 G = origR >> G_SHIFT;
800 B = origR >> B_SHIFT;
801 cachep = &histogram[R * MR + G * MG + B];
802 if (*cachep == 0)
803 fill_inverse_cmap_rgb(quantobj, histogram, R, G, B);
804
805 dest[idx] = quantobj->cmap[*cachep - 1].r;
806
807 /* If the colormap entry for this pixel is the same as the
808 background's colormap entry, set the pixel to the
809 background color. */
810 if (bgColor && dest[idx] == bg_color.r)
811 dest[idx] = bgColor->r;
812 }
813 }
814 }
815
initialize_median_cut(int num_colors)816 static QuantizeObj *initialize_median_cut(int num_colors)
817 {
818 QuantizeObj *quantobj;
819
820 /* Initialize the data structures */
821 XMALLOC(quantobj, sizeof(QuantizeObj));
822
823 XMALLOC(quantobj->histogram, sizeof(ColorFreq) * HIST_R_ELEMS * HIST_G_ELEMS * HIST_B_ELEMS);
824 quantobj->desired_number_of_colors = num_colors;
825
826 return quantobj;
827 }
828
quantize(at_bitmap * image,long ncolors,const at_color * bgColor,QuantizeObj ** iQuant,at_exception_type * exp)829 void quantize(at_bitmap * image, long ncolors, const at_color * bgColor, QuantizeObj ** iQuant, at_exception_type * exp)
830 {
831 QuantizeObj *quantobj;
832 unsigned int spp = AT_BITMAP_PLANES(image);
833
834 if (spp != 3 && spp != 1) {
835 LOG("quantize: %u-plane images are not supported", spp);
836 at_exception_fatal(exp, "quantize: wrong plane images are passed");
837 return;
838 }
839
840 /* If a pointer was sent in, let's use it. */
841 if (iQuant) {
842 if (*iQuant == NULL) {
843 quantobj = initialize_median_cut(ncolors);
844 median_cut_pass1_rgb(quantobj, image, bgColor);
845 *iQuant = quantobj;
846 } else
847 quantobj = *iQuant;
848 } else {
849 quantobj = initialize_median_cut(ncolors);
850 median_cut_pass1_rgb(quantobj, image, NULL);
851 }
852
853 median_cut_pass2_rgb(quantobj, image, bgColor);
854
855 if (iQuant == NULL)
856 quantize_object_free(quantobj);
857 }
858
quantize_object_free(QuantizeObj * quantobj)859 void quantize_object_free(QuantizeObj * quantobj)
860 {
861 free(quantobj->histogram);
862 free(quantobj);
863 }
864