1 : /******************************************************************************
2 : * $Id: gdalmediancut.cpp 13346 2007-12-15 12:53:04Z 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 13346 2007-12-15 12:53:04Z 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 : int CPL_STDCALL
108 : 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 3 : void * pProgressArg )
116 :
117 : {
118 3 : VALIDATE_POINTER1( hRed, "GDALComputeMedianCutPCT", CE_Failure );
119 3 : VALIDATE_POINTER1( hGreen, "GDALComputeMedianCutPCT", CE_Failure );
120 3 : VALIDATE_POINTER1( hBlue, "GDALComputeMedianCutPCT", CE_Failure );
121 :
122 : int nXSize, nYSize;
123 3 : CPLErr err = CE_None;
124 :
125 : /* -------------------------------------------------------------------- */
126 : /* Validate parameters. */
127 : /* -------------------------------------------------------------------- */
128 3 : nXSize = GDALGetRasterBandXSize( hRed );
129 3 : nYSize = GDALGetRasterBandYSize( hRed );
130 :
131 3 : 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 3 : 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 3 : 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 3 : 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 3 : if( pfnProgress == NULL )
168 1 : 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 3 : CPLCalloc(GMC_B_LEN * GMC_B_LEN * GMC_B_LEN,sizeof(int));
181 3 : usedboxes = NULL;
182 3 : box_list = freeboxes = (Colorbox *)CPLMalloc(nColors*sizeof (Colorbox));
183 3 : freeboxes[0].next = &freeboxes[1];
184 3 : freeboxes[0].prev = NULL;
185 277 : for (i = 1; i < nColors-1; ++i) {
186 274 : freeboxes[i].next = &freeboxes[i+1];
187 274 : freeboxes[i].prev = &freeboxes[i-1];
188 : }
189 3 : freeboxes[nColors-1].next = NULL;
190 3 : 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 3 : ptr = freeboxes;
202 3 : freeboxes = ptr->next;
203 3 : if (freeboxes)
204 3 : freeboxes->prev = NULL;
205 3 : ptr->next = usedboxes;
206 3 : usedboxes = ptr;
207 3 : if (ptr->next)
208 0 : ptr->next->prev = ptr;
209 :
210 3 : ptr->rmin = ptr->gmin = ptr->bmin = 999;
211 3 : ptr->rmax = ptr->gmax = ptr->bmax = -1;
212 3 : ptr->total = nXSize * nYSize;
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* Collect histogram. */
216 : /* -------------------------------------------------------------------- */
217 3 : pabyRedLine = (GByte *) VSIMalloc(nXSize);
218 3 : pabyGreenLine = (GByte *) VSIMalloc(nXSize);
219 3 : pabyBlueLine = (GByte *) VSIMalloc(nXSize);
220 :
221 3 : 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 153 : for( iLine = 0; iLine < nYSize; iLine++ )
232 : {
233 150 : 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 150 : pabyRedLine, nXSize, 1, GDT_Byte, 0, 0 );
243 : GDALRasterIO( hGreen, GF_Read, 0, iLine, nXSize, 1,
244 150 : pabyGreenLine, nXSize, 1, GDT_Byte, 0, 0 );
245 : GDALRasterIO( hBlue, GF_Read, 0, iLine, nXSize, 1,
246 150 : pabyBlueLine, nXSize, 1, GDT_Byte, 0, 0 );
247 :
248 7650 : for( iPixel = 0; iPixel < nXSize; iPixel++ )
249 : {
250 : int nRed, nGreen, nBlue;
251 :
252 7500 : nRed = pabyRedLine[iPixel] >> COLOR_SHIFT;
253 7500 : nGreen = pabyGreenLine[iPixel] >> COLOR_SHIFT;
254 7500 : nBlue = pabyBlueLine[iPixel] >> COLOR_SHIFT;
255 :
256 7500 : ptr->rmin = MIN(ptr->rmin, nRed);
257 7500 : ptr->gmin = MIN(ptr->gmin, nGreen);
258 7500 : ptr->bmin = MIN(ptr->bmin, nBlue);
259 7500 : ptr->rmax = MAX(ptr->rmax, nRed);
260 7500 : ptr->gmax = MAX(ptr->gmax, nGreen);
261 7500 : ptr->bmax = MAX(ptr->bmax, nBlue);
262 :
263 7500 : histogram[nRed][nGreen][nBlue]++;
264 : }
265 : }
266 :
267 3 : 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 283 : while (freeboxes != NULL) {
279 277 : ptr = largest_box(usedboxes);
280 277 : if (ptr != NULL)
281 277 : splitbox(ptr, histogram, &freeboxes, &usedboxes);
282 : else
283 0 : freeboxes = NULL;
284 : }
285 :
286 : /* ==================================================================== */
287 : /* STEP 4: assign colors to all boxes */
288 : /* ==================================================================== */
289 283 : for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next)
290 : {
291 : GDALColorEntry sEntry;
292 :
293 280 : sEntry.c1 = (GByte) (((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2);
294 280 : sEntry.c2 = (GByte) (((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2);
295 280 : sEntry.c3 = (GByte) (((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2);
296 280 : sEntry.c4 = 255;
297 280 : GDALSetColorEntry( hColorTable, i, &sEntry );
298 : }
299 :
300 3 : end_and_cleanup:
301 3 : CPLFree( pabyRedLine );
302 3 : CPLFree( pabyGreenLine );
303 3 : CPLFree( pabyBlueLine );
304 :
305 : /* We're done with the boxes now */
306 3 : CPLFree(box_list);
307 3 : freeboxes = usedboxes = NULL;
308 :
309 3 : CPLFree( histogram );
310 :
311 3 : return err;
312 : }
313 :
314 : /************************************************************************/
315 : /* largest_box() */
316 : /************************************************************************/
317 :
318 : static Colorbox *
319 277 : largest_box(Colorbox *usedboxes)
320 : {
321 : Colorbox *p, *b;
322 : int size;
323 :
324 277 : b = NULL;
325 277 : size = -1;
326 33065 : for (p = usedboxes; p != NULL; p = p->next)
327 32788 : if ((p->rmax > p->rmin || p->gmax > p->gmin ||
328 : p->bmax > p->bmin) && p->total > size)
329 1246 : size = (b = p)->total;
330 277 : return (b);
331 : }
332 :
333 : /************************************************************************/
334 : /* splitbox() */
335 : /************************************************************************/
336 : static void
337 : splitbox(Colorbox* ptr, int (*histogram)[GMC_B_LEN][GMC_B_LEN],
338 277 : Colorbox **pfreeboxes, Colorbox **pusedboxes)
339 : {
340 : int hist2[GMC_B_LEN];
341 277 : 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 277 : i = ptr->rmax - ptr->rmin;
355 375 : if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
356 98 : axis = RED;
357 179 : else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
358 99 : axis = GREEN;
359 : else
360 80 : axis = BLUE;
361 : /* get histogram along longest axis */
362 277 : switch (axis) {
363 : case RED:
364 98 : histp = &hist2[ptr->rmin];
365 521 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
366 423 : *histp = 0;
367 4563 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
368 4140 : iptr = &histogram[ir][ig][ptr->bmin];
369 79620 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
370 75480 : *histp += *iptr++;
371 : }
372 423 : histp++;
373 : }
374 98 : first = ptr->rmin;
375 98 : last = ptr->rmax;
376 98 : break;
377 : case GREEN:
378 99 : histp = &hist2[ptr->gmin];
379 618 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
380 519 : *histp = 0;
381 2980 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
382 2461 : iptr = &histogram[ir][ig][ptr->bmin];
383 17083 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
384 14622 : *histp += *iptr++;
385 : }
386 519 : histp++;
387 : }
388 99 : first = ptr->gmin;
389 99 : last = ptr->gmax;
390 99 : break;
391 : case BLUE:
392 80 : histp = &hist2[ptr->bmin];
393 495 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) {
394 415 : *histp = 0;
395 3615 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
396 3200 : iptr = &histogram[ir][ptr->gmin][ib];
397 43860 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
398 40660 : *histp += *iptr;
399 40660 : iptr += GMC_B_LEN;
400 : }
401 : }
402 415 : histp++;
403 : }
404 80 : first = ptr->bmin;
405 80 : last = ptr->bmax;
406 : break;
407 : }
408 : /* find median point */
409 277 : sum2 = ptr->total / 2;
410 277 : histp = &hist2[first];
411 277 : sum = 0;
412 277 : for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
413 : ;
414 277 : if (i == first)
415 114 : i++;
416 :
417 : /* Create new box, re-allocate points */
418 277 : new_cb = *pfreeboxes;
419 277 : *pfreeboxes = new_cb->next;
420 277 : if (*pfreeboxes)
421 274 : (*pfreeboxes)->prev = NULL;
422 277 : if (*pusedboxes)
423 277 : (*pusedboxes)->prev = new_cb;
424 277 : new_cb->next = *pusedboxes;
425 277 : *pusedboxes = new_cb;
426 :
427 277 : histp = &hist2[first];
428 740 : for (sum1 = 0, j = first; j < i; j++)
429 463 : sum1 += *histp++;
430 1171 : for (sum2 = 0, j = i; j <= last; j++)
431 894 : sum2 += *histp++;
432 277 : new_cb->total = sum1;
433 277 : ptr->total = sum2;
434 :
435 277 : new_cb->rmin = ptr->rmin;
436 277 : new_cb->rmax = ptr->rmax;
437 277 : new_cb->gmin = ptr->gmin;
438 277 : new_cb->gmax = ptr->gmax;
439 277 : new_cb->bmin = ptr->bmin;
440 277 : new_cb->bmax = ptr->bmax;
441 277 : switch (axis) {
442 : case RED:
443 98 : new_cb->rmax = i-1;
444 98 : ptr->rmin = i;
445 98 : break;
446 : case GREEN:
447 99 : new_cb->gmax = i-1;
448 99 : ptr->gmin = i;
449 99 : break;
450 : case BLUE:
451 80 : new_cb->bmax = i-1;
452 80 : ptr->bmin = i;
453 : break;
454 : }
455 277 : shrinkbox(new_cb, histogram);
456 277 : shrinkbox(ptr, histogram);
457 277 : }
458 :
459 : /************************************************************************/
460 : /* shrinkbox() */
461 : /************************************************************************/
462 : static void
463 554 : shrinkbox(Colorbox* box, int (*histogram)[GMC_B_LEN][GMC_B_LEN])
464 : {
465 : int *histp, ir, ig, ib;
466 :
467 554 : if (box->rmax > box->rmin) {
468 333 : for (ir = box->rmin; ir <= box->rmax; ++ir)
469 660 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
470 607 : histp = &histogram[ir][ig][box->bmin];
471 3224 : for (ib = box->bmin; ib <= box->bmax; ++ib)
472 2897 : if (*histp++ != 0) {
473 280 : box->rmin = ir;
474 280 : goto have_rmin;
475 : }
476 : }
477 280 : have_rmin:
478 280 : if (box->rmax > box->rmin)
479 450 : for (ir = box->rmax; ir >= box->rmin; --ir)
480 2682 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
481 2503 : histp = &histogram[ir][ig][box->bmin];
482 2503 : ib = box->bmin;
483 14935 : for (; ib <= box->bmax; ++ib)
484 12703 : if (*histp++ != 0) {
485 271 : box->rmax = ir;
486 271 : goto have_rmax;
487 : }
488 : }
489 : }
490 554 : have_rmax:
491 554 : if (box->gmax > box->gmin) {
492 350 : for (ig = box->gmin; ig <= box->gmax; ++ig)
493 1272 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
494 1185 : histp = &histogram[ir][ig][box->bmin];
495 18783 : for (ib = box->bmin; ib <= box->bmax; ++ib)
496 17861 : if (*histp++ != 0) {
497 263 : box->gmin = ig;
498 263 : goto have_gmin;
499 : }
500 : }
501 263 : have_gmin:
502 263 : if (box->gmax > box->gmin)
503 388 : for (ig = box->gmax; ig >= box->gmin; --ig)
504 1437 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
505 1306 : histp = &histogram[ir][ig][box->bmin];
506 1306 : ib = box->bmin;
507 15653 : for (; ib <= box->bmax; ++ib)
508 14604 : if (*histp++ != 0) {
509 257 : box->gmax = ig;
510 257 : goto have_gmax;
511 : }
512 : }
513 : }
514 554 : have_gmax:
515 554 : if (box->bmax > box->bmin) {
516 325 : for (ib = box->bmin; ib <= box->bmax; ++ib)
517 495 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
518 453 : histp = &histogram[ir][box->gmin][ib];
519 1868 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
520 1698 : if (*histp != 0) {
521 283 : box->bmin = ib;
522 283 : goto have_bmin;
523 : }
524 1415 : histp += GMC_B_LEN;
525 : }
526 : }
527 283 : have_bmin:
528 283 : if (box->bmax > box->bmin)
529 432 : for (ib = box->bmax; ib >= box->bmin; --ib)
530 1595 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
531 1427 : histp = &histogram[ir][box->gmin][ib];
532 1427 : ig = box->gmin;
533 16104 : for (; ig <= box->gmax; ++ig) {
534 14941 : if (*histp != 0) {
535 264 : box->bmax = ib;
536 264 : goto have_bmax;
537 : }
538 14677 : histp += GMC_B_LEN;
539 : }
540 : }
541 : }
542 554 : have_bmax:
543 : ;
544 554 : }
|
