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