1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 - This software is distributed in the hope that it will be
4 - useful, but with NO WARRANTY OF ANY KIND.
5 - No author or distributor accepts responsibility to anyone for the
6 - consequences of using this software, or for whether it serves any
7 - particular purpose or works at all, unless he or she says so in
8 - writing. Everyone is granted permission to copy, modify and
9 - redistribute this source code, for commercial or non-commercial
10 - purposes, with the following restrictions: (1) the origin of this
11 - source code must not be misrepresented; (2) modified versions must
12 - be plainly marked as such; and (3) this notice may not be removed
13 - or altered from any source or modified source distribution.
14 *====================================================================*/
15
16 /*
17 * colorquant2.c
18 *
19 * Modified median cut color quantization
20 *
21 * PIX *pixMedianCutQuant()
22 * PIX *pixMedianCutQuantGeneral()
23 * l_int32 *pixMedianCutHisto()
24 *
25 * Static helpers
26 * static PIXCMAP *pixcmapGenerateFromHisto()
27 * static PIX *pixQuantizeWithColormap()
28 * static void getColorIndexMedianCut()
29 * static L_BOX3D *pixGetColorRegion()
30 * static l_int32 medianCutApply()
31 * static PIXCMAP *pixcmapGenerateFromMedianCuts()
32 * static l_int32 vboxGetAverageColor()
33 * static l_int32 vboxGetCount()
34 * static l_int32 vboxGetVolume()
35 * static L_BOX3D *box3dCreate();
36 * static L_BOX3D *box3dCopy();
37 *
38 * Paul Heckbert published the median cut algorithm, "Color Image
39 * Quantization for Frame Buffer Display," in Proc. SIGGRAPH '82,
40 * Boston, July 1982, pp. 297-307. A copy of the paper without
41 * figures can be found on the web.
42 *
43 * Median cut starts with either the full color space or the occupied
44 * region of color space. If you're not dithering, the occupied region
45 * can be used, but with dithering, pixels can end up in any place
46 * in the color space, so you must represent the entire color space in
47 * the final colormap.
48 *
49 * Color components are quantized to typically 5 or 6 significant
50 * bits (for each of r, g and b). Call a 3D region of color
51 * space a 'vbox'. Any color in this quantized space is represented
52 * by an element of a linear histogram array, indexed by rgb value.
53 * The initial region is then divided into two regions that have roughly
54 * equal pixel occupancy (hence the name "median cut"). Subdivision
55 * continues until the requisite number of vboxes has been generated.
56 *
57 * But the devil is in the details of the subdivision process.
58 * Here are some choices that you must make:
59 * (1) Along which axis to subdivide?
60 * (2) Which box to put the bin with the median pixel?
61 * (3) How to order the boxes for subdivision?
62 * (4) How to adequately handle boxes with very small numbers of pixels?
63 * (5) How to prevent a little-represented but highly visible color
64 * from being masked out by other colors in its vbox.
65 *
66 * Taking these in order:
67 * (1) Heckbert suggests using either the largest vbox side, or the vbox
68 * side with the largest variance in pixel occupancy. We choose
69 * to divide based on the largest vbox side.
70 * (2) Suppose you've chosen a side. Then you have a histogram
71 * of pixel occupancy in 2D slices of the vbox. One of those
72 * slices includes the median pixel. Suppose there are L bins
73 * to the left (smaller index) and R bins to the right. Then
74 * this slice (or bin) should be assigned to the box containing
75 * the smaller of L and R. This both shortens the larger
76 * of the subdivided dimensions and helps a low-count color
77 * far from the subdivision boundary to better express itself.
78 * (2a) One can also ask if the boundary should be moved even
79 * farther into the longer side. This is feasable if we have
80 * a method for doing extra subdivisions on the high count
81 * vboxes. And we do (see (3)).
82 * (3) To make sure that the boxes are subdivided toward equal
83 * occupancy, use an occupancy-sorted priority queue, rather
84 * than a simple queue.
85 * (4) With a priority queue, boxes with small number of pixels
86 * won't be repeatedly subdivided. This is good.
87 * (5) Use of a priority queue allows tricks such as in (2a) to let
88 * small occupancy clusters be better expressed. In addition,
89 * rather than splitting near the median, small occupancy colors
90 * are best reproduced by cutting half-way into the longer side.
91 *
92 * However, serious problems can arise with dithering if a priority
93 * queue is used based on population alone. If the picture has
94 * large regions of nearly constant color, some vboxes can be very
95 * large and have a sizeable population (but not big enough to get to
96 * the head of the queue). If one of these large, occupied vboxes
97 * is near in color to a nearly constant color region of the
98 * image, dithering can inject pixels from the large vbox into
99 * the nearly uniform region. These pixels can be very far away
100 * in color, and the oscillations are highly visible. To prevent
101 * this, we can take either or both of these actions:
102 *
103 * (1) Subdivide a fraction (< 1.0) based on population, and
104 * do the rest of the subdivision based on the product of
105 * the vbox volume and its population. By using the product,
106 * we avoid further subdivision of nearly empty vboxes, and
107 * directly target large vboxes with significant population.
108 *
109 * (2) Threshold the excess color transferred in dithering to
110 * neighboring pixels.
111 *
112 * Doing either of these will stop the most annoying oscillations
113 * in dithering. Furthermore, by doing (1), we also improve the
114 * rendering of regions of nearly constant color, both with and
115 * without dithering. It turns out that the image quality is
116 * not sensitive to the value of the parameter in (1); values
117 * between 0.3 and 0.9 give very good results.
118 *
119 * Here's the lesson: subdivide the color space into vboxes such
120 * that (1) the most populated vboxes that can be further
121 * subdivided (i.e., that occupy more than one quantum volume
122 * in color space) all have approximately the same population,
123 * and (2) all large vboxes have no significant population.
124 * If these conditions are met, the quantization will be excellent.
125 *
126 * Once the subdivision has been made, the colormap is generated,
127 * with one color for each vbox and using the average color in the vbox.
128 * At the same time, the histogram array is converted to an inverse
129 * colormap table, storing the colormap index in every cell in the
130 * vbox. Finally, using both the colormap and the inverse colormap,
131 * a colormapped pix is quickly generated from the original rgb pix.
132 *
133 * In the present implementation, subdivided regions of colorspace
134 * that are not occupied are retained, but not further subdivided.
135 * This is required for our inverse colormap lookup table for
136 * dithering, because dithered pixels may fall into these unoccupied
137 * regions. For such empty regions, we use the center as the rgb
138 * colormap value.
139 *
140 * This variation on median cut can be referred to as "Modified Median
141 * Cut" quantization, or MMCQ. Overall, the undithered MMCQ gives
142 * comparable results to the two-pass Octcube Quantizer (OQ).
143 * Comparing the two methods on the test24.jpg painting, we see:
144 *
145 * (1) For rendering spot color (the various reds and pinks in
146 * the image), MMCQ is not as good as OQ.
147 *
148 * (2) For rendering majority color regions, MMCQ does a better
149 * job of avoiding posterization. That is, it does better
150 * dividing the color space up in the most heavily populated regions.
151 */
152
153 #include <stdio.h>
154 #include <stdlib.h>
155 #include <string.h>
156 #include <math.h>
157 #include "allheaders.h"
158
159 /* Median cut 3-d volume element. Sort on first element, which
160 * can be the number of pixels, the volume or a combination
161 * of these. */
162 struct L_Box3d
163 {
164 l_float32 sortparam; /* parameter on which to sort the vbox */
165 l_int32 npix; /* number of pixels in the vbox */
166 l_int32 vol; /* quantized volume of vbox */
167 l_int32 r1; /* min r index in the vbox */
168 l_int32 r2; /* max r index in the vbox */
169 l_int32 g1; /* min g index in the vbox */
170 l_int32 g2; /* max g index in the vbox */
171 l_int32 b1; /* min b index in the vbox */
172 l_int32 b2; /* max b index in the vbox */
173 };
174 typedef struct L_Box3d L_BOX3D;
175
176 /* Static median cut helper functions */
177 static PIXCMAP *pixcmapGenerateFromHisto(PIX *pixs, l_int32 depth,
178 l_int32 *histo, l_int32 histosize,
179 l_int32 sigbits);
180 static PIX *pixQuantizeWithColormap(PIX *pixs, l_int32 ditherflag,
181 l_int32 outdepth,
182 PIXCMAP *cmap, l_int32 *indexmap,
183 l_int32 mapsize, l_int32 sigbits);
184 static void getColorIndexMedianCut(l_uint32 pixel, l_int32 rshift,
185 l_uint32 mask, l_int32 sigbits,
186 l_int32 *pindex);
187 static L_BOX3D *pixGetColorRegion(PIX *pixs, l_int32 sigbits,
188 l_int32 subsample);
189 static l_int32 medianCutApply(l_int32 *histo, l_int32 sigbits,
190 L_BOX3D *vbox, L_BOX3D **pvbox1,
191 L_BOX3D **pvbox2);
192 static PIXCMAP *pixcmapGenerateFromMedianCuts(L_HEAP *lh, l_int32 *histo,
193 l_int32 sigbits);
194 static l_int32 vboxGetAverageColor(L_BOX3D *vbox, l_int32 *histo,
195 l_int32 sigbits, l_int32 index,
196 l_int32 *prval, l_int32 *pgval,
197 l_int32 *pbval);
198 static l_int32 vboxGetCount(L_BOX3D *vbox, l_int32 *histo, l_int32 sigbits);
199 static l_int32 vboxGetVolume(L_BOX3D *vbox);
200 static L_BOX3D *box3dCreate(l_int32 r1, l_int32 r2, l_int32 g1,
201 l_int32 g2, l_int32 b1, l_int32 b2);
202 static L_BOX3D *box3dCopy(L_BOX3D *vbox);
203
204
205 /* 5 significant bits for each component is generally satisfactory */
206 static const l_int32 DEFAULT_SIG_BITS = 5;
207 static const l_int32 MAX_ITERS_ALLOWED = 5000; /* prevents infinite looping */
208
209 /* Specify fraction of vboxes made that are sorted on population alone.
210 * The remaining vboxes are sorted on (population * vbox-volume). */
211 static const l_float32 FRACT_BY_POPULATION = 0.85;
212
213 /* To get the max value of 'dif' in the dithering color transfer,
214 * divide DIF_CAP by 8. */
215 static const l_int32 DIF_CAP = 100;
216
217
218 #ifndef NO_CONSOLE_IO
219 #define DEBUG_MC_COLORS 0
220 #define DEBUG_SPLIT_AXES 0
221 #endif /* ~NO_CONSOLE_IO */
222
223
224 /*!
225 * pixMedianCutQuant()
226 *
227 * Input: pixs (32 bpp; rgb color)
228 * ditherflag (1 for dither; 0 for no dither)
229 * Return: pixd (8 bit with colormap), or null on error
230 *
231 * Notes:
232 * (1) Simple interface. See pixMedianCutQuantGeneral() for
233 * use of defaulted parameters.
234 */
235 PIX *
pixMedianCutQuant(PIX * pixs,l_int32 ditherflag)236 pixMedianCutQuant(PIX *pixs,
237 l_int32 ditherflag)
238 {
239 return pixMedianCutQuantGeneral(pixs, ditherflag,
240 0, 256, DEFAULT_SIG_BITS, 1);
241 }
242
243
244 /*!
245 * pixMedianCutQuantGeneral()
246 *
247 * Input: pixs (32 bpp; rgb color)
248 * ditherflag (1 for dither; 0 for no dither)
249 * outdepth (output depth; valid: 0, 1, 2, 4, 8)
250 * maxcolors (between 2 and 256)
251 * sigbits (valid: 5 or 6; use 0 for default)
252 * maxsub (max subsampling, integer; use 0 for default;
253 * 1 for no subsampling)
254 * Return: pixd (8 bit with colormap), or null on error
255 *
256 * Notes:
257 * (1) @maxcolors must be in the range [2 ... 256].
258 * (2) Use @outdepth = 0 to have the output depth computed as the
259 * minimum required to hold the actual colors found, given
260 * the @maxcolors constraint.
261 * (3) Use @outdepth = 1, 2, 4 or 8 to specify the output depth.
262 * In that case, @maxcolors must not exceed 2^(outdepth).
263 * (4) If there are fewer quantized colors in the image than @maxcolors,
264 * the colormap is simply generated from those colors.
265 * (5) @maxsub is the maximum allowed subsampling to be used in the
266 * computation of the color histogram and region of occupied
267 * color space. The subsampling is chosen internally for
268 * efficiency, based on the image size, but this parameter
269 * limits it. Use @maxsub = 0 for the internal default, which is the
270 * maximum allowed subsampling. Use @maxsub = 1 to prevent
271 * subsampling. In general use @maxsub >= 1 to specify the
272 * maximum subsampling to be allowed, where the actual subsampling
273 * will be the minimum of this value and the internally
274 * determined default value.
275 * (6) If the image appears gray because either most of the pixels
276 * are gray or most of the pixels are essentially black or white,
277 * the image is trivially quantized with a grayscale colormap. The
278 * reason is that median cut divides the color space into rectangular
279 * regions, and it does a very poor job if all the pixels are
280 * near the diagonal of the color space cube.
281 */
282 PIX *
pixMedianCutQuantGeneral(PIX * pixs,l_int32 ditherflag,l_int32 outdepth,l_int32 maxcolors,l_int32 sigbits,l_int32 maxsub)283 pixMedianCutQuantGeneral(PIX *pixs,
284 l_int32 ditherflag,
285 l_int32 outdepth,
286 l_int32 maxcolors,
287 l_int32 sigbits,
288 l_int32 maxsub)
289 {
290 l_int32 i, subsample, histosize, smalln, ncolors, niters, popcolors;
291 l_int32 w, h, minside, factor;
292 l_int32 *histo;
293 l_float32 pixfract, colorfract;
294 L_BOX3D *vbox, *vbox1, *vbox2;
295 L_HEAP *lh, *lhs;
296 PIX *pixd;
297 PIXCMAP *cmap;
298
299 PROCNAME("pixMedianCutQuantGeneral");
300
301 if (!pixs || pixGetDepth(pixs) != 32)
302 return (PIX *)ERROR_PTR("pixs undefined or not 32 bpp", procName, NULL);
303 if (maxcolors < 2 || maxcolors > 256)
304 return (PIX *)ERROR_PTR("maxcolors not in [2...256]", procName, NULL);
305 if (outdepth != 0 && outdepth != 1 && outdepth != 2 && outdepth != 4 &&
306 outdepth != 8)
307 return (PIX *)ERROR_PTR("outdepth not in {0,1,2,4,8}", procName, NULL);
308 if (outdepth > 0 && (maxcolors > (1 << outdepth)))
309 return (PIX *)ERROR_PTR("maxcolors > 2^(outdepth)", procName, NULL);
310 if (sigbits == 0)
311 sigbits = DEFAULT_SIG_BITS;
312 else if (sigbits < 5 || sigbits > 6)
313 return (PIX *)ERROR_PTR("sigbits not 5 or 6", procName, NULL);
314 if (maxsub <= 0)
315 maxsub = 10; /* default will prevail for 10^7 pixels or less */
316
317 /* Determine if the image has sufficient color content.
318 * If pixfract << 1, most pixels are close to black or white.
319 * If colorfract << 1, the pixels that are not near
320 * black or white have very little color.
321 * If without color, quantize with a grayscale colormap. */
322 pixGetDimensions(pixs, &w, &h, NULL);
323 minside = L_MIN(w, h);
324 factor = L_MAX(1, minside / 200);
325 pixColorFraction(pixs, 20, 248, 12, factor, &pixfract, &colorfract);
326 if (pixfract < 0.03 || colorfract < 0.03) {
327 L_INFO_FLOAT2("\n Pixel fraction neither white nor black = %6.3f"
328 "\n Color fraction of those pixels = %6.3f"
329 "\n Quantizing in gray",
330 procName, pixfract, colorfract);
331 return pixConvertTo8(pixs, 1);
332 }
333
334 /* Compute the color space histogram. Default sampling
335 * is about 10^5 pixels. */
336 if (maxsub == 1)
337 subsample = 1;
338 else {
339 subsample = (l_int32)(sqrt((l_float64)(w * h) / 100000.));
340 subsample = L_MAX(1, L_MIN(maxsub, subsample));
341 }
342 histo = pixMedianCutHisto(pixs, sigbits, subsample);
343 histosize = 1 << (3 * sigbits);
344
345 /* See if the number of quantized colors is less than maxcolors */
346 ncolors = 0;
347 smalln = TRUE;
348 for (i = 0; i < histosize; i++) {
349 if (histo[i])
350 ncolors++;
351 if (ncolors > maxcolors) {
352 smalln = FALSE;
353 break;
354 }
355 }
356 if (smalln) { /* finish up now */
357 if (outdepth == 0) {
358 if (ncolors <= 2)
359 outdepth = 1;
360 else if (ncolors <= 4)
361 outdepth = 2;
362 else if (ncolors <= 16)
363 outdepth = 4;
364 else
365 outdepth = 8;
366 }
367 cmap = pixcmapGenerateFromHisto(pixs, outdepth,
368 histo, histosize, sigbits);
369 pixd = pixQuantizeWithColormap(pixs, ditherflag, outdepth, cmap,
370 histo, histosize, sigbits);
371 FREE(histo);
372 return pixd;
373 }
374
375 /* Initial vbox: minimum region in colorspace occupied by pixels */
376 if (ditherflag || subsample > 1) /* use full color space */
377 vbox = box3dCreate(0, (1 << sigbits) - 1,
378 0, (1 << sigbits) - 1,
379 0, (1 << sigbits) - 1);
380 else
381 vbox = pixGetColorRegion(pixs, sigbits, subsample);
382 vbox->npix = vboxGetCount(vbox, histo, sigbits);
383 vbox->vol = vboxGetVolume(vbox);
384
385 /* For a fraction 'popcolors' of the desired 'maxcolors',
386 * generate median cuts based on population, putting
387 * everything on a priority queue sorted by population. */
388 lh = lheapCreate(0, L_SORT_DECREASING);
389 lheapAdd(lh, vbox);
390 ncolors = 1;
391 niters = 0;
392 popcolors = (l_int32)(FRACT_BY_POPULATION * maxcolors);
393 while (1) {
394 vbox = (L_BOX3D *)lheapRemove(lh);
395 if (vboxGetCount(vbox, histo, sigbits) == 0) { /* just put it back */
396 lheapAdd(lh, vbox);
397 continue;
398 }
399 medianCutApply(histo, sigbits, vbox, &vbox1, &vbox2);
400 if (!vbox1) {
401 L_WARNING("vbox1 not defined; shouldn't happen!", procName);
402 break;
403 }
404 if (vbox1->vol > 1)
405 vbox1->sortparam = vbox1->npix;
406 FREE(vbox);
407 lheapAdd(lh, vbox1);
408 if (vbox2) { /* vbox2 can be NULL */
409 if (vbox2->vol > 1)
410 vbox2->sortparam = vbox2->npix;
411 lheapAdd(lh, vbox2);
412 ncolors++;
413 }
414 if (ncolors >= popcolors)
415 break;
416 if (niters++ > MAX_ITERS_ALLOWED) {
417 L_WARNING("infinite loop; perhaps too few pixels!", procName);
418 break;
419 }
420 }
421
422 /* Re-sort by the product of pixel occupancy times the size
423 * in color space. */
424 lhs = lheapCreate(0, L_SORT_DECREASING);
425 while ((vbox = (L_BOX3D *)lheapRemove(lh))) {
426 vbox->sortparam = vbox->npix * vbox->vol;
427 lheapAdd(lhs, vbox);
428 }
429 lheapDestroy(&lh, TRUE);
430
431 /* For the remaining (maxcolors - popcolors), generate the
432 * median cuts using the (npix * vol) sorting. */
433 while (1) {
434 vbox = (L_BOX3D *)lheapRemove(lhs);
435 if (vboxGetCount(vbox, histo, sigbits) == 0) { /* just put it back */
436 lheapAdd(lhs, vbox);
437 continue;
438 }
439 medianCutApply(histo, sigbits, vbox, &vbox1, &vbox2);
440 if (!vbox1) {
441 L_WARNING("vbox1 not defined; shouldn't happen!", procName);
442 break;
443 }
444 if (vbox1->vol > 1)
445 vbox1->sortparam = vbox1->npix * vbox1->vol;
446 FREE(vbox);
447 lheapAdd(lhs, vbox1);
448 if (vbox2) { /* vbox2 can be NULL */
449 if (vbox2->vol > 1)
450 vbox2->sortparam = vbox2->npix * vbox2->vol;
451 lheapAdd(lhs, vbox2);
452 ncolors++;
453 }
454 if (ncolors >= maxcolors)
455 break;
456 if (niters++ > MAX_ITERS_ALLOWED) {
457 L_WARNING("infinite loop; perhaps too few pixels!", procName);
458 break;
459 }
460 }
461
462 /* Re-sort by pixel occupancy. This is not necessary,
463 * but it makes a more useful listing. */
464 lh = lheapCreate(0, L_SORT_DECREASING);
465 while ((vbox = (L_BOX3D *)lheapRemove(lhs))) {
466 vbox->sortparam = vbox->npix;
467 /* vbox->sortparam = vbox->npix * vbox->vol; */
468 lheapAdd(lh, vbox);
469 }
470 lheapDestroy(&lhs, TRUE);
471
472 /* Generate colormap from median cuts and quantize pixd */
473 cmap = pixcmapGenerateFromMedianCuts(lh, histo, sigbits);
474 if (outdepth == 0) {
475 ncolors = pixcmapGetCount(cmap);
476 if (ncolors <= 2)
477 outdepth = 1;
478 else if (ncolors <= 4)
479 outdepth = 2;
480 else if (ncolors <= 16)
481 outdepth = 4;
482 else
483 outdepth = 8;
484 }
485 pixd = pixQuantizeWithColormap(pixs, ditherflag, outdepth, cmap,
486 histo, histosize, sigbits);
487
488 lheapDestroy(&lh, TRUE);
489 FREE(histo);
490 return pixd;
491 }
492
493
494 /*!
495 * pixMedianCutHisto()
496 *
497 * Input: pixs (32 bpp; rgb color)
498 * sigbits (valid: 5 or 6)
499 * subsample (integer > 0)
500 * Return: histo (1-d array, giving the number of pixels in
501 * each quantized region of color space), or null on error
502 *
503 * Notes:
504 * (1) Array is indexed by (3 * sigbits) bits. The array size
505 * is 2^(3 * sigbits).
506 * (2) Indexing into the array from rgb uses red sigbits as
507 * most significant and blue as least.
508 */
509 l_int32 *
pixMedianCutHisto(PIX * pixs,l_int32 sigbits,l_int32 subsample)510 pixMedianCutHisto(PIX *pixs,
511 l_int32 sigbits,
512 l_int32 subsample)
513 {
514 l_int32 i, j, w, h, wpl, rshift, index, histosize;
515 l_int32 *histo;
516 l_uint32 mask, pixel;
517 l_uint32 *data, *line;
518
519 PROCNAME("pixMedianCutHisto");
520
521 if (!pixs)
522 return (l_int32 *)ERROR_PTR("pixs not defined", procName, NULL);
523 if (pixGetDepth(pixs) != 32)
524 return (l_int32 *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
525 if (sigbits < 5 || sigbits > 6)
526 return (l_int32 *)ERROR_PTR("sigbits not 5 or 6", procName, NULL);
527 if (subsample <= 0)
528 return (l_int32 *)ERROR_PTR("subsample not > 0", procName, NULL);
529
530 histosize = 1 << (3 * sigbits);
531 if ((histo = (l_int32 *)CALLOC(histosize, sizeof(l_int32))) == NULL)
532 return (l_int32 *)ERROR_PTR("histo not made", procName, NULL);
533
534 rshift = 8 - sigbits;
535 mask = 0xff >> rshift;
536 pixGetDimensions(pixs, &w, &h, NULL);
537 data = pixGetData(pixs);
538 wpl = pixGetWpl(pixs);
539 for (i = 0; i < h; i += subsample) {
540 line = data + i * wpl;
541 for (j = 0; j < w; j += subsample) {
542 pixel = line[j];
543 getColorIndexMedianCut(pixel, rshift, mask, sigbits, &index);
544 histo[index]++;
545 }
546 }
547
548 return histo;
549 }
550
551
552 /*!
553 * pixcmapGenerateFromHisto()
554 *
555 * Input: pixs (32 bpp; rgb color)
556 * depth (of colormap)
557 * histo
558 * histosize
559 * sigbits
560 * Return: colormap, or null on error
561 *
562 * Notes:
563 * (1) This is used when the number of colors in the histo
564 * is not greater than maxcolors.
565 * (2) As a side-effect, the histo becomes an inverse colormap,
566 * labeling the cmap indices for each existing color.
567 */
568 static PIXCMAP *
pixcmapGenerateFromHisto(PIX * pixs,l_int32 depth,l_int32 * histo,l_int32 histosize,l_int32 sigbits)569 pixcmapGenerateFromHisto(PIX *pixs,
570 l_int32 depth,
571 l_int32 *histo,
572 l_int32 histosize,
573 l_int32 sigbits)
574 {
575 l_int32 i, index, shift, rval, gval, bval;
576 l_uint32 mask;
577 PIXCMAP *cmap;
578
579 PROCNAME("pixcmapGenerateFromHisto");
580
581 if (!pixs)
582 return (PIXCMAP *)ERROR_PTR("pixs not defined", procName, NULL);
583 if (pixGetDepth(pixs) != 32)
584 return (PIXCMAP *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
585 if (!histo)
586 return (PIXCMAP *)ERROR_PTR("histo not defined", procName, NULL);
587
588 /* Capture the rgb values of each occupied cube in the histo,
589 * and re-label the histo value with the colormap index. */
590 cmap = pixcmapCreate(depth);
591 shift = 8 - sigbits;
592 mask = 0xff >> shift;
593 for (i = 0, index = 0; i < histosize; i++) {
594 if (histo[i]) {
595 rval = (i >> (2 * sigbits)) << shift;
596 gval = ((i >> sigbits) & mask) << shift;
597 bval = (i & mask) << shift;
598 pixcmapAddColor(cmap, rval, gval, bval);
599 histo[i] = index++;
600 }
601 }
602
603 return cmap;
604 }
605
606
607 /*!
608 * pixQuantizeWithColormap()
609 *
610 * Input: pixs (32 bpp; rgb color)
611 * ditherflag (1 for dither; 0 for no dither)
612 * outdepth
613 * cmap
614 * indexmap
615 * histosize
616 * sigbits
617 * Return: pixd (quantized to colormap), or null on error
618 *
619 * Notes:
620 * (1) The indexmap is a LUT that takes the rgb indices of the
621 * pixel and returns the index into the colormap.
622 * (2) If ditherflag is 1, @outdepth is ignored and the output
623 * depth is set to 8.
624 */
625 static PIX *
pixQuantizeWithColormap(PIX * pixs,l_int32 ditherflag,l_int32 outdepth,PIXCMAP * cmap,l_int32 * indexmap,l_int32 mapsize,l_int32 sigbits)626 pixQuantizeWithColormap(PIX *pixs,
627 l_int32 ditherflag,
628 l_int32 outdepth,
629 PIXCMAP *cmap,
630 l_int32 *indexmap,
631 l_int32 mapsize,
632 l_int32 sigbits)
633 {
634 l_uint8 *bufu8r, *bufu8g, *bufu8b;
635 l_int32 i, j, w, h, wpls, wpld, rshift, index, cmapindex;
636 l_int32 rval, gval, bval, rc, gc, bc;
637 l_int32 dif, val1, val2, val3;
638 l_int32 *buf1r, *buf1g, *buf1b, *buf2r, *buf2g, *buf2b;
639 l_uint32 *datas, *datad, *lines, *lined;
640 l_uint32 mask, pixel;
641 PIX *pixd;
642
643 PROCNAME("pixQuantizeWithColormap");
644
645 if (!pixs || pixGetDepth(pixs) != 32)
646 return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL);
647 if (!cmap)
648 return (PIX *)ERROR_PTR("cmap not defined", procName, NULL);
649 if (!indexmap)
650 return (PIX *)ERROR_PTR("indexmap not defined", procName, NULL);
651 if (ditherflag)
652 outdepth = 8;
653
654 pixGetDimensions(pixs, &w, &h, NULL);
655 pixd = pixCreate(w, h, outdepth);
656 pixSetColormap(pixd, cmap);
657 pixCopyResolution(pixd, pixs);
658 pixCopyInputFormat(pixd, pixs);
659 datas = pixGetData(pixs);
660 datad = pixGetData(pixd);
661 wpls = pixGetWpl(pixs);
662 wpld = pixGetWpl(pixd);
663
664 rshift = 8 - sigbits;
665 mask = 0xff >> rshift;
666 if (ditherflag == 0) {
667 for (i = 0; i < h; i++) {
668 lines = datas + i * wpls;
669 lined = datad + i * wpld;
670 if (outdepth == 1) {
671 for (j = 0; j < w; j++) {
672 pixel = lines[j];
673 getColorIndexMedianCut(pixel, rshift, mask,
674 sigbits, &index);
675 if (indexmap[index])
676 SET_DATA_BIT(lined, j);
677 }
678 }
679 else if (outdepth == 2) {
680 for (j = 0; j < w; j++) {
681 pixel = lines[j];
682 getColorIndexMedianCut(pixel, rshift, mask,
683 sigbits, &index);
684 SET_DATA_DIBIT(lined, j, indexmap[index]);
685 }
686 }
687 else if (outdepth == 4) {
688 for (j = 0; j < w; j++) {
689 pixel = lines[j];
690 getColorIndexMedianCut(pixel, rshift, mask,
691 sigbits, &index);
692 SET_DATA_QBIT(lined, j, indexmap[index]);
693 }
694 }
695 else { /* outdepth == 8 */
696 for (j = 0; j < w; j++) {
697 pixel = lines[j];
698 getColorIndexMedianCut(pixel, rshift, mask,
699 sigbits, &index);
700 SET_DATA_BYTE(lined, j, indexmap[index]);
701 }
702 }
703 }
704 }
705 else { /* ditherflag == 1 */
706 bufu8r = (l_uint8 *)CALLOC(w, sizeof(l_uint8));
707 bufu8g = (l_uint8 *)CALLOC(w, sizeof(l_uint8));
708 bufu8b = (l_uint8 *)CALLOC(w, sizeof(l_uint8));
709 buf1r = (l_int32 *)CALLOC(w, sizeof(l_int32));
710 buf1g = (l_int32 *)CALLOC(w, sizeof(l_int32));
711 buf1b = (l_int32 *)CALLOC(w, sizeof(l_int32));
712 buf2r = (l_int32 *)CALLOC(w, sizeof(l_int32));
713 buf2g = (l_int32 *)CALLOC(w, sizeof(l_int32));
714 buf2b = (l_int32 *)CALLOC(w, sizeof(l_int32));
715 if (!bufu8r || !bufu8g || !bufu8b)
716 return (PIX *)ERROR_PTR("uint8 line buf not made", procName, NULL);
717 if (!buf1r || !buf1g || !buf1b || !buf2r || !buf2g || !buf2b)
718 return (PIX *)ERROR_PTR("mono line buf not made", procName, NULL);
719
720 /* Start by priming buf2; line 1 is above line 2 */
721 pixGetRGBLine(pixs, 0, bufu8r, bufu8g, bufu8b);
722 for (j = 0; j < w; j++) {
723 buf2r[j] = 64 * bufu8r[j];
724 buf2g[j] = 64 * bufu8g[j];
725 buf2b[j] = 64 * bufu8b[j];
726 }
727
728 for (i = 0; i < h - 1; i++) {
729 /* Swap data 2 --> 1, and read in new line 2 */
730 memcpy(buf1r, buf2r, 4 * w);
731 memcpy(buf1g, buf2g, 4 * w);
732 memcpy(buf1b, buf2b, 4 * w);
733 pixGetRGBLine(pixs, i + 1, bufu8r, bufu8g, bufu8b);
734 for (j = 0; j < w; j++) {
735 buf2r[j] = 64 * bufu8r[j];
736 buf2g[j] = 64 * bufu8g[j];
737 buf2b[j] = 64 * bufu8b[j];
738 }
739
740 /* Dither */
741 lined = datad + i * wpld;
742 for (j = 0; j < w - 1; j++) {
743 rval = buf1r[j] / 64;
744 gval = buf1g[j] / 64;
745 bval = buf1b[j] / 64;
746 index = ((rval >> rshift) << (2 * sigbits)) +
747 ((gval >> rshift) << sigbits) + (bval >> rshift);
748 cmapindex = indexmap[index];
749 SET_DATA_BYTE(lined, j, cmapindex);
750 pixcmapGetColor(cmap, cmapindex, &rc, &gc, &bc);
751
752 dif = buf1r[j] / 8 - 8 * rc;
753 if (dif > DIF_CAP) dif = DIF_CAP;
754 if (dif < -DIF_CAP) dif = -DIF_CAP;
755 if (dif != 0) {
756 val1 = buf1r[j + 1] + 3 * dif;
757 val2 = buf2r[j] + 3 * dif;
758 val3 = buf2r[j + 1] + 2 * dif;
759 if (dif > 0) {
760 buf1r[j + 1] = L_MIN(16383, val1);
761 buf2r[j] = L_MIN(16383, val2);
762 buf2r[j + 1] = L_MIN(16383, val3);
763 }
764 else if (dif < 0) {
765 buf1r[j + 1] = L_MAX(0, val1);
766 buf2r[j] = L_MAX(0, val2);
767 buf2r[j + 1] = L_MAX(0, val3);
768 }
769 }
770
771 dif = buf1g[j] / 8 - 8 * gc;
772 if (dif > DIF_CAP) dif = DIF_CAP;
773 if (dif < -DIF_CAP) dif = -DIF_CAP;
774 if (dif != 0) {
775 val1 = buf1g[j + 1] + 3 * dif;
776 val2 = buf2g[j] + 3 * dif;
777 val3 = buf2g[j + 1] + 2 * dif;
778 if (dif > 0) {
779 buf1g[j + 1] = L_MIN(16383, val1);
780 buf2g[j] = L_MIN(16383, val2);
781 buf2g[j + 1] = L_MIN(16383, val3);
782 }
783 else if (dif < 0) {
784 buf1g[j + 1] = L_MAX(0, val1);
785 buf2g[j] = L_MAX(0, val2);
786 buf2g[j + 1] = L_MAX(0, val3);
787 }
788 }
789
790 dif = buf1b[j] / 8 - 8 * bc;
791 if (dif > DIF_CAP) dif = DIF_CAP;
792 if (dif < -DIF_CAP) dif = -DIF_CAP;
793 if (dif != 0) {
794 val1 = buf1b[j + 1] + 3 * dif;
795 val2 = buf2b[j] + 3 * dif;
796 val3 = buf2b[j + 1] + 2 * dif;
797 if (dif > 0) {
798 buf1b[j + 1] = L_MIN(16383, val1);
799 buf2b[j] = L_MIN(16383, val2);
800 buf2b[j + 1] = L_MIN(16383, val3);
801 }
802 else if (dif < 0) {
803 buf1b[j + 1] = L_MAX(0, val1);
804 buf2b[j] = L_MAX(0, val2);
805 buf2b[j + 1] = L_MAX(0, val3);
806 }
807 }
808 }
809
810 /* Get last pixel in row; no downward propagation */
811 rval = buf1r[w - 1] / 64;
812 gval = buf1g[w - 1] / 64;
813 bval = buf1b[w - 1] / 64;
814 index = ((rval >> rshift) << (2 * sigbits)) +
815 ((gval >> rshift) << sigbits) + (bval >> rshift);
816 SET_DATA_BYTE(lined, w - 1, indexmap[index]);
817 }
818
819 /* Get last row of pixels; no leftward propagation */
820 lined = datad + (h - 1) * wpld;
821 for (j = 0; j < w; j++) {
822 rval = buf2r[j] / 64;
823 gval = buf2g[j] / 64;
824 bval = buf2b[j] / 64;
825 index = ((rval >> rshift) << (2 * sigbits)) +
826 ((gval >> rshift) << sigbits) + (bval >> rshift);
827 SET_DATA_BYTE(lined, j, indexmap[index]);
828 }
829
830 FREE(bufu8r);
831 FREE(bufu8g);
832 FREE(bufu8b);
833 FREE(buf1r);
834 FREE(buf1g);
835 FREE(buf1b);
836 FREE(buf2r);
837 FREE(buf2g);
838 FREE(buf2b);
839 }
840
841 return pixd;
842 }
843
844
845 /*!
846 * getColorIndexMedianCut()
847 *
848 * Input: pixel (32 bit rgb)
849 * rshift (of component: 8 - sigbits)
850 * mask (over sigbits)
851 * sigbits
852 * &index (<return> rgb index value)
853 * Return: void
854 *
855 * Notes:
856 * (1) This is used on each pixel in the source image. No checking
857 * is done on input values.
858 */
859 static void
getColorIndexMedianCut(l_uint32 pixel,l_int32 rshift,l_uint32 mask,l_int32 sigbits,l_int32 * pindex)860 getColorIndexMedianCut(l_uint32 pixel,
861 l_int32 rshift,
862 l_uint32 mask,
863 l_int32 sigbits,
864 l_int32 *pindex)
865 {
866 l_int32 rval, gval, bval;
867
868 rval = pixel >> (24 + rshift);
869 gval = (pixel >> (16 + rshift)) & mask;
870 bval = (pixel >> (8 + rshift)) & mask;
871 *pindex = (rval << (2 * sigbits)) + (gval << sigbits) + bval;
872 return;
873 }
874
875
876 /*!
877 * pixGetColorRegion()
878 *
879 * Input: pixs (32 bpp; rgb color)
880 * sigbits (valid: 5, 6)
881 * subsample (integer > 0)
882 * Return: vbox (minimum 3D box in color space enclosing all pixels),
883 * or null on error
884 *
885 * Notes:
886 * (1) Computes the minimum 3D box in color space enclosing all
887 * pixels in the image.
888 */
889 static L_BOX3D *
pixGetColorRegion(PIX * pixs,l_int32 sigbits,l_int32 subsample)890 pixGetColorRegion(PIX *pixs,
891 l_int32 sigbits,
892 l_int32 subsample)
893 {
894 l_int32 rmin, rmax, gmin, gmax, bmin, bmax, rval, gval, bval;
895 l_int32 w, h, wpl, i, j, rshift;
896 l_uint32 mask, pixel;
897 l_uint32 *data, *line;
898
899 PROCNAME("pixGetColorRegion");
900
901 if (!pixs)
902 return (L_BOX3D *)ERROR_PTR("pixs not defined", procName, NULL);
903
904 rmin = gmin = bmin = 1000000;
905 rmax = gmax = bmax = 0;
906 rshift = 8 - sigbits;
907 mask = 0xff >> rshift;
908 pixGetDimensions(pixs, &w, &h, NULL);
909 data = pixGetData(pixs);
910 wpl = pixGetWpl(pixs);
911 for (i = 0; i < h; i += subsample) {
912 line = data + i * wpl;
913 for (j = 0; j < w; j += subsample) {
914 pixel = line[j];
915 rval = pixel >> (24 + rshift);
916 gval = (pixel >> (16 + rshift)) & mask;
917 bval = (pixel >> (8 + rshift)) & mask;
918 if (rval < rmin)
919 rmin = rval;
920 else if (rval > rmax)
921 rmax = rval;
922 if (gval < gmin)
923 gmin = gval;
924 else if (gval > gmax)
925 gmax = gval;
926 if (bval < bmin)
927 bmin = bval;
928 else if (bval > bmax)
929 bmax = bval;
930 }
931 }
932
933 return box3dCreate(rmin, rmax, gmin, gmax, bmin, bmax);
934 }
935
936
937 /*!
938 * medianCutApply()
939 *
940 * Input: histo (array; in rgb colorspace)
941 * sigbits
942 * vbox (input 3D box)
943 * &vbox1, vbox2 (<return> vbox split in two parts)
944 * Return: 0 if OK, 1 on error
945 */
946 static l_int32
medianCutApply(l_int32 * histo,l_int32 sigbits,L_BOX3D * vbox,L_BOX3D ** pvbox1,L_BOX3D ** pvbox2)947 medianCutApply(l_int32 *histo,
948 l_int32 sigbits,
949 L_BOX3D *vbox,
950 L_BOX3D **pvbox1,
951 L_BOX3D **pvbox2)
952 {
953 l_int32 i, j, k, sum, rw, gw, bw, maxw, index;
954 l_int32 total, left, right;
955 l_int32 partialsum[128];
956 L_BOX3D *vbox1, *vbox2;
957
958 PROCNAME("medianCutApply");
959
960 if (!histo)
961 return ERROR_INT("histo not defined", procName, 1);
962 if (!vbox)
963 return ERROR_INT("vbox not defined", procName, 1);
964 if (!pvbox1 || !pvbox2)
965 return ERROR_INT("&vbox1 and &vbox2 not both defined", procName, 1);
966
967 *pvbox1 = *pvbox2 = NULL;
968 if (vboxGetCount(vbox, histo, sigbits) == 0)
969 return ERROR_INT("no pixels in vbox", procName, 1);
970
971 /* If the vbox occupies just one element in color space, it can't
972 * be split. Leave the 'sortparam' field at 0, so that it goes to
973 * the tail of the priority queue and stays there, thereby avoiding
974 * an infinite loop (take off, put back on the head) if it
975 * happens to be the most populous box! */
976 rw = vbox->r2 - vbox->r1 + 1;
977 gw = vbox->g2 - vbox->g1 + 1;
978 bw = vbox->b2 - vbox->b1 + 1;
979 if (rw == 1 && gw == 1 && bw == 1) {
980 *pvbox1 = box3dCopy(vbox);
981 return 0;
982 }
983
984 /* Select the longest axis for splitting */
985 maxw = L_MAX(rw, gw);
986 maxw = L_MAX(maxw, bw);
987 #if DEBUG_SPLIT_AXES
988 if (rw == maxw)
989 fprintf(stderr, "red split\n");
990 else if (gw == maxw)
991 fprintf(stderr, "green split\n");
992 else
993 fprintf(stderr, "blue split\n");
994 #endif /* DEBUG_SPLIT_AXES */
995
996 /* Find the partial sum arrays along the selected axis. */
997 total = 0;
998 if (maxw == rw) {
999 for (i = vbox->r1; i <= vbox->r2; i++) {
1000 sum = 0;
1001 for (j = vbox->g1; j <= vbox->g2; j++) {
1002 for (k = vbox->b1; k <= vbox->b2; k++) {
1003 index = (i << (2 * sigbits)) + (j << sigbits) + k;
1004 sum += histo[index];
1005 }
1006 }
1007 total += sum;
1008 partialsum[i] = total;
1009 }
1010 }
1011 else if (maxw == gw) {
1012 for (i = vbox->g1; i <= vbox->g2; i++) {
1013 sum = 0;
1014 for (j = vbox->r1; j <= vbox->r2; j++) {
1015 for (k = vbox->b1; k <= vbox->b2; k++) {
1016 index = (i << sigbits) + (j << (2 * sigbits)) + k;
1017 sum += histo[index];
1018 }
1019 }
1020 total += sum;
1021 partialsum[i] = total;
1022 }
1023 }
1024 else { /* maxw == bw */
1025 for (i = vbox->b1; i <= vbox->b2; i++) {
1026 sum = 0;
1027 for (j = vbox->r1; j <= vbox->r2; j++) {
1028 for (k = vbox->g1; k <= vbox->g2; k++) {
1029 index = i + (j << (2 * sigbits)) + (k << sigbits);
1030 sum += histo[index];
1031 }
1032 }
1033 total += sum;
1034 partialsum[i] = total;
1035 }
1036 }
1037
1038 /* Determine the cut planes, making sure that two vboxes
1039 * are always produced. Generate the two vboxes and compute
1040 * the sum in each of them. Choose the cut plane within
1041 * the greater of the (left, right) sides of the bin in which
1042 * the median pixel resides. Here's the surprise: go halfway
1043 * into that side. By doing that, you technically move away
1044 * from "median cut," but in the process a significant number
1045 * of low-count vboxes are produced, allowing much better
1046 * reproduction of low-count spot colors. */
1047 if (maxw == rw) {
1048 for (i = vbox->r1; i <= vbox->r2; i++) {
1049 if (partialsum[i] > total / 2) {
1050 vbox1 = box3dCopy(vbox);
1051 vbox2 = box3dCopy(vbox);
1052 left = i - vbox->r1;
1053 right = vbox->r2 - i;
1054 if (left <= right)
1055 vbox1->r2 = L_MIN(vbox->r2 - 1, i + right / 2);
1056 else /* left > right */
1057 vbox1->r2 = L_MAX(vbox->r1, i - 1 - left / 2);
1058 vbox2->r1 = vbox1->r2 + 1;
1059 break;
1060 }
1061 }
1062 }
1063 else if (maxw == gw) {
1064 for (i = vbox->g1; i <= vbox->g2; i++) {
1065 if (partialsum[i] > total / 2) {
1066 vbox1 = box3dCopy(vbox);
1067 vbox2 = box3dCopy(vbox);
1068 left = i - vbox->g1;
1069 right = vbox->g2 - i;
1070 if (left <= right)
1071 vbox1->g2 = L_MIN(vbox->g2 - 1, i + right / 2);
1072 else /* left > right */
1073 vbox1->g2 = L_MAX(vbox->g1, i - 1 - left / 2);
1074 vbox2->g1 = vbox1->g2 + 1;
1075 break;
1076 }
1077 }
1078 }
1079 else { /* maxw == bw */
1080 for (i = vbox->b1; i <= vbox->b2; i++) {
1081 if (partialsum[i] > total / 2) {
1082 vbox1 = box3dCopy(vbox);
1083 vbox2 = box3dCopy(vbox);
1084 left = i - vbox->b1;
1085 right = vbox->b2 - i;
1086 if (left <= right)
1087 vbox1->b2 = L_MIN(vbox->b2 - 1, i + right / 2);
1088 else /* left > right */
1089 vbox1->b2 = L_MAX(vbox->b1, i - 1 - left / 2);
1090 vbox2->b1 = vbox1->b2 + 1;
1091 break;
1092 }
1093 }
1094 }
1095 vbox1->npix = vboxGetCount(vbox1, histo, sigbits);
1096 vbox2->npix = vboxGetCount(vbox2, histo, sigbits);
1097 vbox1->vol = vboxGetVolume(vbox1);
1098 vbox2->vol = vboxGetVolume(vbox2);
1099 *pvbox1 = vbox1;
1100 *pvbox2 = vbox2;
1101
1102 return 0;
1103 }
1104
1105
1106 /*!
1107 * pixcmapGenerateFromMedianCuts()
1108 *
1109 * Input: lh (priority queue of pointers to vboxes)
1110 * histo
1111 * sigbits (valid: 5 or 6)
1112 * Return: cmap, or null on error
1113 *
1114 * Notes:
1115 * (1) Each vbox in the heap represents a color in the colormap.
1116 * (2) As a side-effect, the histo becomes an inverse colormap,
1117 * where the part of the array correpsonding to each vbox
1118 * is labeled with the cmap index for that vbox. Then
1119 * for each rgb pixel, the colormap index is found directly
1120 * by mapping the rgb value to the histo array index.
1121 */
1122 static PIXCMAP *
pixcmapGenerateFromMedianCuts(L_HEAP * lh,l_int32 * histo,l_int32 sigbits)1123 pixcmapGenerateFromMedianCuts(L_HEAP *lh,
1124 l_int32 *histo,
1125 l_int32 sigbits)
1126 {
1127 l_int32 index, rval, gval, bval;
1128 L_BOX3D *vbox;
1129 PIXCMAP *cmap;
1130
1131 PROCNAME("pixcmapGenerateFromMedianCuts");
1132
1133 if (!lh)
1134 return (PIXCMAP *)ERROR_PTR("lh not defined", procName, NULL);
1135 if (!histo)
1136 return (PIXCMAP *)ERROR_PTR("histo not defined", procName, NULL);
1137
1138 rval = gval = bval = 0; /* make compiler happy */
1139 cmap = pixcmapCreate(8);
1140 index = 0;
1141 while (lheapGetCount(lh) > 0) {
1142 vbox = (L_BOX3D *)lheapRemove(lh);
1143 vboxGetAverageColor(vbox, histo, sigbits, index, &rval, &gval, &bval);
1144 pixcmapAddColor(cmap, rval, gval, bval);
1145 FREE(vbox);
1146 index++;
1147 }
1148
1149 return cmap;
1150 }
1151
1152
1153 /*!
1154 * vboxGetAverageColor()
1155 *
1156 * Input: vbox (3d region of color space for one quantized color)
1157 * histo
1158 * sigbits (valid: 5 or 6)
1159 * index (if >= 0, assign to all colors in histo in this vbox)
1160 * &rval, &gval, &bval (<returned> average color)
1161 * Return: cmap, or null on error
1162 *
1163 * Notes:
1164 * (1) The vbox represents one color in the colormap.
1165 * (2) If index >= 0, as a side-effect, all array elements in
1166 * the histo corresponding to the vbox are labeled with this
1167 * cmap index for that vbox. Otherwise, the histo array
1168 * is not changed.
1169 */
1170 static l_int32
vboxGetAverageColor(L_BOX3D * vbox,l_int32 * histo,l_int32 sigbits,l_int32 index,l_int32 * prval,l_int32 * pgval,l_int32 * pbval)1171 vboxGetAverageColor(L_BOX3D *vbox,
1172 l_int32 *histo,
1173 l_int32 sigbits,
1174 l_int32 index,
1175 l_int32 *prval,
1176 l_int32 *pgval,
1177 l_int32 *pbval)
1178 {
1179 l_int32 i, j, k, ntot, mult, histoindex, rsum, gsum, bsum;
1180
1181 PROCNAME("vboxGetAverageColor");
1182
1183 if (!vbox)
1184 return ERROR_INT("vbox not defined", procName, 1);
1185 if (!histo)
1186 return ERROR_INT("histo not defined", procName, 1);
1187 if (!prval || !pgval || !pbval)
1188 return ERROR_INT("&p*val not all defined", procName, 1);
1189
1190 *prval = *pgval = *pbval = 0;
1191 ntot = 0;
1192 mult = 1 << (8 - sigbits);
1193 rsum = gsum = bsum = 0;
1194 for (i = vbox->r1; i <= vbox->r2; i++) {
1195 for (j = vbox->g1; j <= vbox->g2; j++) {
1196 for (k = vbox->b1; k <= vbox->b2; k++) {
1197 histoindex = (i << (2 * sigbits)) + (j << sigbits) + k;
1198 ntot += histo[histoindex];
1199 rsum += (l_int32)(histo[histoindex] * (i + 0.5) * mult);
1200 gsum += (l_int32)(histo[histoindex] * (j + 0.5) * mult);
1201 bsum += (l_int32)(histo[histoindex] * (k + 0.5) * mult);
1202 if (index >= 0)
1203 histo[histoindex] = index;
1204 }
1205 }
1206 }
1207
1208 if (ntot == 0) {
1209 *prval = mult * (vbox->r1 + vbox->r2 + 1) / 2;
1210 *pgval = mult * (vbox->g1 + vbox->g2 + 1) / 2;
1211 *pbval = mult * (vbox->b1 + vbox->b2 + 1) / 2;
1212 }
1213 else {
1214 *prval = rsum / ntot;
1215 *pgval = gsum / ntot;
1216 *pbval = bsum / ntot;
1217 }
1218
1219 #if DEBUG_MC_COLORS
1220 fprintf(stderr, "ntot[%d] = %d: [%d, %d, %d], (%d, %d, %d)\n",
1221 index, ntot, vbox->r2 - vbox->r1 + 1,
1222 vbox->g2 - vbox->g1 + 1, vbox->b2 - vbox->b1 + 1,
1223 *prval, *pgval, *pbval);
1224 #endif /* DEBUG_MC_COLORS */
1225
1226 return 0;
1227 }
1228
1229
1230 /*!
1231 * vboxGetCount()
1232 *
1233 * Input: vbox (3d region of color space for one quantized color)
1234 * histo
1235 * sigbits (valid: 5 or 6)
1236 * Return: number of image pixels in this region, or 0 on error
1237 */
1238 static l_int32
vboxGetCount(L_BOX3D * vbox,l_int32 * histo,l_int32 sigbits)1239 vboxGetCount(L_BOX3D *vbox,
1240 l_int32 *histo,
1241 l_int32 sigbits)
1242 {
1243 l_int32 i, j, k, npix, index;
1244
1245 PROCNAME("vboxGetCount");
1246
1247 if (!vbox)
1248 return ERROR_INT("vbox not defined", procName, 0);
1249 if (!histo)
1250 return ERROR_INT("histo not defined", procName, 0);
1251
1252 npix = 0;
1253 for (i = vbox->r1; i <= vbox->r2; i++) {
1254 for (j = vbox->g1; j <= vbox->g2; j++) {
1255 for (k = vbox->b1; k <= vbox->b2; k++) {
1256 index = (i << (2 * sigbits)) + (j << sigbits) + k;
1257 npix += histo[index];
1258 }
1259 }
1260 }
1261
1262 return npix;
1263 }
1264
1265
1266 /*!
1267 * vboxGetVolume()
1268 *
1269 * Input: vbox (3d region of color space for one quantized color)
1270 * Return: quantized volume of vbox, or 0 on error
1271 */
1272 static l_int32
vboxGetVolume(L_BOX3D * vbox)1273 vboxGetVolume(L_BOX3D *vbox)
1274 {
1275 PROCNAME("vboxGetVolume");
1276
1277 if (!vbox)
1278 return ERROR_INT("vbox not defined", procName, 0);
1279
1280 return ((vbox->r2 - vbox->r1 + 1) * (vbox->g2 - vbox->g1 + 1) *
1281 (vbox->b2 - vbox->b1 + 1));
1282 }
1283
1284
1285 static L_BOX3D *
box3dCreate(l_int32 r1,l_int32 r2,l_int32 g1,l_int32 g2,l_int32 b1,l_int32 b2)1286 box3dCreate(l_int32 r1,
1287 l_int32 r2,
1288 l_int32 g1,
1289 l_int32 g2,
1290 l_int32 b1,
1291 l_int32 b2)
1292 {
1293 L_BOX3D *vbox;
1294
1295 vbox = (L_BOX3D *)CALLOC(1, sizeof(L_BOX3D));
1296 vbox->r1 = r1;
1297 vbox->r2 = r2;
1298 vbox->g1 = g1;
1299 vbox->g2 = g2;
1300 vbox->b1 = b1;
1301 vbox->b2 = b2;
1302 return vbox;
1303 }
1304
1305
1306 /*
1307 * Note: don't copy the sortparam.
1308 */
1309 static L_BOX3D *
box3dCopy(L_BOX3D * vbox)1310 box3dCopy(L_BOX3D *vbox)
1311 {
1312 L_BOX3D *vboxc;
1313
1314 PROCNAME("box3dCopy");
1315
1316 if (!vbox)
1317 return (L_BOX3D *)ERROR_PTR("vbox not defined", procName, NULL);
1318
1319 vboxc = box3dCreate(vbox->r1, vbox->r2, vbox->g1, vbox->g2,
1320 vbox->b1, vbox->b2);
1321 vboxc->npix = vbox->npix;
1322 vboxc->vol = vbox->vol;
1323 return vboxc;
1324 }
1325
1326
1327