1 : /******************************************************************************
2 : * $Id: gdalmediancut.cpp 20629 2010-09-16 19:08:40Z rouault $
3 : *
4 : * Project: CIETMap Phase 2
5 : * Purpose: Use median cut algorithm to generate an near-optimal PCT for a
6 : * given RGB image. Implemented as function GDALComputeMedianCutPCT.
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2001, Frank Warmerdam
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ******************************************************************************
30 : *
31 : * This code was based on the tiffmedian.c code from libtiff (www.libtiff.org)
32 : * which was based on a paper by Paul Heckbert:
33 : *
34 : * "Color Image Quantization for Frame Buffer Display", Paul
35 : * Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307.
36 : *
37 : */
38 :
39 : #include "gdal_priv.h"
40 : #include "gdal_alg.h"
41 :
42 : CPL_CVSID("$Id: gdalmediancut.cpp 20629 2010-09-16 19:08:40Z rouault $");
43 :
44 : #define MAX_CMAP_SIZE 256
45 :
46 : #define COLOR_DEPTH 8
47 : #define MAX_COLOR 256
48 :
49 : #define GMC_B_DEPTH 5 /* # bits/pixel to use */
50 : #define GMC_B_LEN (1L<<GMC_B_DEPTH)
51 :
52 : #define COLOR_SHIFT (COLOR_DEPTH-GMC_B_DEPTH)
53 :
54 : typedef struct colorbox {
55 : struct colorbox *next, *prev;
56 : int rmin, rmax;
57 : int gmin, gmax;
58 : int bmin, bmax;
59 : int total;
60 : } Colorbox;
61 :
62 : static void splitbox(Colorbox* ptr, int (*histogram)[GMC_B_LEN][GMC_B_LEN],
63 : Colorbox **pfreeboxes, Colorbox **pusedboxes);
64 : static void shrinkbox(Colorbox* box, int (*histogram)[GMC_B_LEN][GMC_B_LEN]);
65 : static Colorbox* largest_box(Colorbox *usedboxes);
66 :
67 : /************************************************************************/
68 : /* GDALComputeMedianCutPCT() */
69 : /************************************************************************/
70 :
71 : /**
72 : * Compute optimal PCT for RGB image.
73 : *
74 : * This function implements a median cut algorithm to compute an "optimal"
75 : * pseudocolor table for representing an input RGB image. This PCT could
76 : * then be used with GDALDitherRGB2PCT() to convert a 24bit RGB image into
77 : * an eightbit pseudo-colored image.
78 : *
79 : * This code was based on the tiffmedian.c code from libtiff (www.libtiff.org)
80 : * which was based on a paper by Paul Heckbert:
81 : *
82 : * \verbatim
83 : * "Color Image Quantization for Frame Buffer Display", Paul
84 : * Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307.
85 : * \endverbatim
86 : *
87 : * The red, green and blue input bands do not necessarily need to come
88 : * from the same file, but they must be the same width and height. They will
89 : * be clipped to 8bit during reading, so non-eight bit bands are generally
90 : * inappropriate.
91 : *
92 : * @param hRed Red input band.
93 : * @param hGreen Green input band.
94 : * @param hBlue Blue input band.
95 : * @param pfnIncludePixel function used to test which pixels should be included
96 : * in the analysis. At this time this argument is ignored and all pixels are
97 : * utilized. This should normally be NULL.
98 : * @param nColors the desired number of colors to be returned (2-256).
99 : * @param hColorTable the colors will be returned in this color table object.
100 : * @param pfnProgress callback for reporting algorithm progress matching the
101 : * GDALProgressFunc() semantics. May be NULL.
102 : * @param pProgressArg callback argument passed to pfnProgress.
103 : *
104 : * @return returns CE_None on success or CE_Failure if an error occurs.
105 : */
106 :
107 : extern "C" int CPL_STDCALL
108 6 : GDALComputeMedianCutPCT( GDALRasterBandH hRed,
109 : GDALRasterBandH hGreen,
110 : GDALRasterBandH hBlue,
111 : int (*pfnIncludePixel)(int,int,void*),
112 : int nColors,
113 : GDALColorTableH hColorTable,
114 : GDALProgressFunc pfnProgress,
115 : void * pProgressArg )
116 :
117 : {
118 6 : VALIDATE_POINTER1( hRed, "GDALComputeMedianCutPCT", CE_Failure );
119 6 : VALIDATE_POINTER1( hGreen, "GDALComputeMedianCutPCT", CE_Failure );
120 6 : VALIDATE_POINTER1( hBlue, "GDALComputeMedianCutPCT", CE_Failure );
121 :
122 : int nXSize, nYSize;
123 6 : CPLErr err = CE_None;
124 :
125 : /* -------------------------------------------------------------------- */
126 : /* Validate parameters. */
127 : /* -------------------------------------------------------------------- */
128 6 : nXSize = GDALGetRasterBandXSize( hRed );
129 6 : nYSize = GDALGetRasterBandYSize( hRed );
130 :
131 6 : if( GDALGetRasterBandXSize( hGreen ) != nXSize
132 : || GDALGetRasterBandYSize( hGreen ) != nYSize
133 : || GDALGetRasterBandXSize( hBlue ) != nXSize
134 : || GDALGetRasterBandYSize( hBlue ) != nYSize )
135 : {
136 : CPLError( CE_Failure, CPLE_IllegalArg,
137 0 : "Green or blue band doesn't match size of red band.\n" );
138 :
139 0 : return CE_Failure;
140 : }
141 :
142 6 : if( pfnIncludePixel != NULL )
143 : {
144 : CPLError( CE_Failure, CPLE_IllegalArg,
145 : "GDALComputeMedianCutPCT() doesn't currently support "
146 0 : " pfnIncludePixel function." );
147 :
148 0 : return CE_Failure;
149 : }
150 :
151 6 : if ( nColors <= 0 )
152 : {
153 : CPLError( CE_Failure, CPLE_IllegalArg,
154 0 : "GDALComputeMedianCutPCT() : nColors must be strictly greater than 1." );
155 :
156 0 : return CE_Failure;
157 : }
158 :
159 6 : if ( nColors > 256 )
160 : {
161 : CPLError( CE_Failure, CPLE_IllegalArg,
162 0 : "GDALComputeMedianCutPCT() : nColors must be lesser than or equal to 256." );
163 :
164 0 : return CE_Failure;
165 : }
166 :
167 6 : if( pfnProgress == NULL )
168 2 : pfnProgress = GDALDummyProgress;
169 :
170 : /* ==================================================================== */
171 : /* STEP 1: crate empty boxes. */
172 : /* ==================================================================== */
173 : int i;
174 : Colorbox *box_list, *ptr;
175 : int (*histogram)[GMC_B_LEN][GMC_B_LEN];
176 : Colorbox *freeboxes;
177 : Colorbox *usedboxes;
178 :
179 : histogram = (int (*)[GMC_B_LEN][GMC_B_LEN])
180 6 : CPLCalloc(GMC_B_LEN * GMC_B_LEN * GMC_B_LEN,sizeof(int));
181 6 : usedboxes = NULL;
182 6 : box_list = freeboxes = (Colorbox *)CPLMalloc(nColors*sizeof (Colorbox));
183 6 : freeboxes[0].next = &freeboxes[1];
184 6 : freeboxes[0].prev = NULL;
185 554 : for (i = 1; i < nColors-1; ++i) {
186 548 : freeboxes[i].next = &freeboxes[i+1];
187 548 : freeboxes[i].prev = &freeboxes[i-1];
188 : }
189 6 : freeboxes[nColors-1].next = NULL;
190 6 : freeboxes[nColors-1].prev = &freeboxes[nColors-2];
191 :
192 : /* ==================================================================== */
193 : /* Build histogram. */
194 : /* ==================================================================== */
195 : GByte *pabyRedLine, *pabyGreenLine, *pabyBlueLine;
196 : int iLine, iPixel;
197 :
198 : /* -------------------------------------------------------------------- */
199 : /* Initialize the box datastructures. */
200 : /* -------------------------------------------------------------------- */
201 6 : ptr = freeboxes;
202 6 : freeboxes = ptr->next;
203 6 : if (freeboxes)
204 6 : freeboxes->prev = NULL;
205 6 : ptr->next = usedboxes;
206 6 : usedboxes = ptr;
207 6 : if (ptr->next)
208 0 : ptr->next->prev = ptr;
209 :
210 6 : ptr->rmin = ptr->gmin = ptr->bmin = 999;
211 6 : ptr->rmax = ptr->gmax = ptr->bmax = -1;
212 6 : ptr->total = nXSize * nYSize;
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* Collect histogram. */
216 : /* -------------------------------------------------------------------- */
217 6 : pabyRedLine = (GByte *) VSIMalloc(nXSize);
218 6 : pabyGreenLine = (GByte *) VSIMalloc(nXSize);
219 6 : pabyBlueLine = (GByte *) VSIMalloc(nXSize);
220 :
221 6 : if (pabyRedLine == NULL ||
222 : pabyGreenLine == NULL ||
223 : pabyBlueLine == NULL)
224 : {
225 : CPLError( CE_Failure, CPLE_OutOfMemory,
226 0 : "VSIMalloc(): Out of memory in GDALComputeMedianCutPCT" );
227 0 : err = CE_Failure;
228 0 : goto end_and_cleanup;
229 : }
230 :
231 306 : for( iLine = 0; iLine < nYSize; iLine++ )
232 : {
233 300 : if( !pfnProgress( iLine / (double) nYSize,
234 : "Generating Histogram", pProgressArg ) )
235 : {
236 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
237 0 : err = CE_Failure;
238 0 : goto end_and_cleanup;
239 : }
240 :
241 : GDALRasterIO( hRed, GF_Read, 0, iLine, nXSize, 1,
242 300 : pabyRedLine, nXSize, 1, GDT_Byte, 0, 0 );
243 : GDALRasterIO( hGreen, GF_Read, 0, iLine, nXSize, 1,
244 300 : pabyGreenLine, nXSize, 1, GDT_Byte, 0, 0 );
245 : GDALRasterIO( hBlue, GF_Read, 0, iLine, nXSize, 1,
246 300 : pabyBlueLine, nXSize, 1, GDT_Byte, 0, 0 );
247 :
248 15300 : for( iPixel = 0; iPixel < nXSize; iPixel++ )
249 : {
250 : int nRed, nGreen, nBlue;
251 :
252 15000 : nRed = pabyRedLine[iPixel] >> COLOR_SHIFT;
253 15000 : nGreen = pabyGreenLine[iPixel] >> COLOR_SHIFT;
254 15000 : nBlue = pabyBlueLine[iPixel] >> COLOR_SHIFT;
255 :
256 15000 : ptr->rmin = MIN(ptr->rmin, nRed);
257 15000 : ptr->gmin = MIN(ptr->gmin, nGreen);
258 15000 : ptr->bmin = MIN(ptr->bmin, nBlue);
259 15000 : ptr->rmax = MAX(ptr->rmax, nRed);
260 15000 : ptr->gmax = MAX(ptr->gmax, nGreen);
261 15000 : ptr->bmax = MAX(ptr->bmax, nBlue);
262 :
263 15000 : histogram[nRed][nGreen][nBlue]++;
264 : }
265 : }
266 :
267 6 : if( !pfnProgress( 1.0, "Generating Histogram", pProgressArg ) )
268 : {
269 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
270 0 : err = CE_Failure;
271 0 : goto end_and_cleanup;
272 : }
273 :
274 : /* ==================================================================== */
275 : /* STEP 3: continually subdivide boxes until no more free */
276 : /* boxes remain or until all colors assigned. */
277 : /* ==================================================================== */
278 566 : while (freeboxes != NULL) {
279 554 : ptr = largest_box(usedboxes);
280 554 : if (ptr != NULL)
281 554 : splitbox(ptr, histogram, &freeboxes, &usedboxes);
282 : else
283 0 : freeboxes = NULL;
284 : }
285 :
286 : /* ==================================================================== */
287 : /* STEP 4: assign colors to all boxes */
288 : /* ==================================================================== */
289 566 : for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next)
290 : {
291 : GDALColorEntry sEntry;
292 :
293 560 : sEntry.c1 = (GByte) (((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2);
294 560 : sEntry.c2 = (GByte) (((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2);
295 560 : sEntry.c3 = (GByte) (((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2);
296 560 : sEntry.c4 = 255;
297 560 : GDALSetColorEntry( hColorTable, i, &sEntry );
298 : }
299 :
300 : end_and_cleanup:
301 6 : CPLFree( pabyRedLine );
302 6 : CPLFree( pabyGreenLine );
303 6 : CPLFree( pabyBlueLine );
304 :
305 : /* We're done with the boxes now */
306 6 : CPLFree(box_list);
307 6 : freeboxes = usedboxes = NULL;
308 :
309 6 : CPLFree( histogram );
310 :
311 6 : return err;
312 : }
313 :
314 : /************************************************************************/
315 : /* largest_box() */
316 : /************************************************************************/
317 :
318 : static Colorbox *
319 554 : largest_box(Colorbox *usedboxes)
320 : {
321 : Colorbox *p, *b;
322 : int size;
323 :
324 554 : b = NULL;
325 554 : size = -1;
326 66130 : for (p = usedboxes; p != NULL; p = p->next)
327 65576 : if ((p->rmax > p->rmin || p->gmax > p->gmin ||
328 : p->bmax > p->bmin) && p->total > size)
329 2492 : size = (b = p)->total;
330 554 : return (b);
331 : }
332 :
333 : /************************************************************************/
334 : /* splitbox() */
335 : /************************************************************************/
336 : static void
337 554 : splitbox(Colorbox* ptr, int (*histogram)[GMC_B_LEN][GMC_B_LEN],
338 : Colorbox **pfreeboxes, Colorbox **pusedboxes)
339 : {
340 : int hist2[GMC_B_LEN];
341 554 : int first=0, last=0;
342 : Colorbox *new_cb;
343 : int *iptr, *histp;
344 : int i, j;
345 : int ir,ig,ib;
346 : int sum, sum1, sum2;
347 : enum { RED, GREEN, BLUE } axis;
348 :
349 : /*
350 : * See which axis is the largest, do a histogram along that
351 : * axis. Split at median point. Contract both new boxes to
352 : * fit points and return
353 : */
354 554 : i = ptr->rmax - ptr->rmin;
355 750 : if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
356 196 : axis = RED;
357 358 : else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
358 198 : axis = GREEN;
359 : else
360 160 : axis = BLUE;
361 : /* get histogram along longest axis */
362 554 : switch (axis) {
363 : case RED:
364 196 : histp = &hist2[ptr->rmin];
365 1042 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
366 846 : *histp = 0;
367 9126 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
368 8280 : iptr = &histogram[ir][ig][ptr->bmin];
369 159240 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
370 150960 : *histp += *iptr++;
371 : }
372 846 : histp++;
373 : }
374 196 : first = ptr->rmin;
375 196 : last = ptr->rmax;
376 196 : break;
377 : case GREEN:
378 198 : histp = &hist2[ptr->gmin];
379 1236 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
380 1038 : *histp = 0;
381 5960 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
382 4922 : iptr = &histogram[ir][ig][ptr->bmin];
383 34166 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
384 29244 : *histp += *iptr++;
385 : }
386 1038 : histp++;
387 : }
388 198 : first = ptr->gmin;
389 198 : last = ptr->gmax;
390 198 : break;
391 : case BLUE:
392 160 : histp = &hist2[ptr->bmin];
393 990 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) {
394 830 : *histp = 0;
395 7230 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
396 6400 : iptr = &histogram[ir][ptr->gmin][ib];
397 87720 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
398 81320 : *histp += *iptr;
399 81320 : iptr += GMC_B_LEN;
400 : }
401 : }
402 830 : histp++;
403 : }
404 160 : first = ptr->bmin;
405 160 : last = ptr->bmax;
406 : break;
407 : }
408 : /* find median point */
409 554 : sum2 = ptr->total / 2;
410 554 : histp = &hist2[first];
411 554 : sum = 0;
412 554 : for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
413 : ;
414 554 : if (i == first)
415 228 : i++;
416 :
417 : /* Create new box, re-allocate points */
418 554 : new_cb = *pfreeboxes;
419 554 : *pfreeboxes = new_cb->next;
420 554 : if (*pfreeboxes)
421 548 : (*pfreeboxes)->prev = NULL;
422 554 : if (*pusedboxes)
423 554 : (*pusedboxes)->prev = new_cb;
424 554 : new_cb->next = *pusedboxes;
425 554 : *pusedboxes = new_cb;
426 :
427 554 : histp = &hist2[first];
428 1480 : for (sum1 = 0, j = first; j < i; j++)
429 926 : sum1 += *histp++;
430 2342 : for (sum2 = 0, j = i; j <= last; j++)
431 1788 : sum2 += *histp++;
432 554 : new_cb->total = sum1;
433 554 : ptr->total = sum2;
434 :
435 554 : new_cb->rmin = ptr->rmin;
436 554 : new_cb->rmax = ptr->rmax;
437 554 : new_cb->gmin = ptr->gmin;
438 554 : new_cb->gmax = ptr->gmax;
439 554 : new_cb->bmin = ptr->bmin;
440 554 : new_cb->bmax = ptr->bmax;
441 554 : switch (axis) {
442 : case RED:
443 196 : new_cb->rmax = i-1;
444 196 : ptr->rmin = i;
445 196 : break;
446 : case GREEN:
447 198 : new_cb->gmax = i-1;
448 198 : ptr->gmin = i;
449 198 : break;
450 : case BLUE:
451 160 : new_cb->bmax = i-1;
452 160 : ptr->bmin = i;
453 : break;
454 : }
455 554 : shrinkbox(new_cb, histogram);
456 554 : shrinkbox(ptr, histogram);
457 554 : }
458 :
459 : /************************************************************************/
460 : /* shrinkbox() */
461 : /************************************************************************/
462 : static void
463 1108 : shrinkbox(Colorbox* box, int (*histogram)[GMC_B_LEN][GMC_B_LEN])
464 : {
465 : int *histp, ir, ig, ib;
466 :
467 1108 : if (box->rmax > box->rmin) {
468 666 : for (ir = box->rmin; ir <= box->rmax; ++ir)
469 1320 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
470 1214 : histp = &histogram[ir][ig][box->bmin];
471 6448 : for (ib = box->bmin; ib <= box->bmax; ++ib)
472 5794 : if (*histp++ != 0) {
473 560 : box->rmin = ir;
474 560 : goto have_rmin;
475 : }
476 : }
477 : have_rmin:
478 560 : if (box->rmax > box->rmin)
479 900 : for (ir = box->rmax; ir >= box->rmin; --ir)
480 5364 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
481 5006 : histp = &histogram[ir][ig][box->bmin];
482 5006 : ib = box->bmin;
483 29870 : for (; ib <= box->bmax; ++ib)
484 25406 : if (*histp++ != 0) {
485 542 : box->rmax = ir;
486 542 : goto have_rmax;
487 : }
488 : }
489 : }
490 : have_rmax:
491 1108 : if (box->gmax > box->gmin) {
492 700 : for (ig = box->gmin; ig <= box->gmax; ++ig)
493 2544 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
494 2370 : histp = &histogram[ir][ig][box->bmin];
495 37566 : for (ib = box->bmin; ib <= box->bmax; ++ib)
496 35722 : if (*histp++ != 0) {
497 526 : box->gmin = ig;
498 526 : goto have_gmin;
499 : }
500 : }
501 : have_gmin:
502 526 : if (box->gmax > box->gmin)
503 776 : for (ig = box->gmax; ig >= box->gmin; --ig)
504 2874 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
505 2612 : histp = &histogram[ir][ig][box->bmin];
506 2612 : ib = box->bmin;
507 31306 : for (; ib <= box->bmax; ++ib)
508 29208 : if (*histp++ != 0) {
509 514 : box->gmax = ig;
510 514 : goto have_gmax;
511 : }
512 : }
513 : }
514 : have_gmax:
515 1108 : if (box->bmax > box->bmin) {
516 650 : for (ib = box->bmin; ib <= box->bmax; ++ib)
517 990 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
518 906 : histp = &histogram[ir][box->gmin][ib];
519 3736 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
520 3396 : if (*histp != 0) {
521 566 : box->bmin = ib;
522 566 : goto have_bmin;
523 : }
524 2830 : histp += GMC_B_LEN;
525 : }
526 : }
527 : have_bmin:
528 566 : if (box->bmax > box->bmin)
529 864 : for (ib = box->bmax; ib >= box->bmin; --ib)
530 3190 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
531 2854 : histp = &histogram[ir][box->gmin][ib];
532 2854 : ig = box->gmin;
533 32208 : for (; ig <= box->gmax; ++ig) {
534 29882 : if (*histp != 0) {
535 528 : box->bmax = ib;
536 528 : goto have_bmax;
537 : }
538 29354 : histp += GMC_B_LEN;
539 : }
540 : }
541 : }
542 : have_bmax:
543 : ;
544 1108 : }
|
