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 1 : 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 1 : VALIDATE_POINTER1( hRed, "GDALComputeMedianCutPCT", CE_Failure );
119 1 : VALIDATE_POINTER1( hGreen, "GDALComputeMedianCutPCT", CE_Failure );
120 1 : VALIDATE_POINTER1( hBlue, "GDALComputeMedianCutPCT", CE_Failure );
121 :
122 : int nXSize, nYSize;
123 1 : CPLErr err = CE_None;
124 :
125 : /* -------------------------------------------------------------------- */
126 : /* Validate parameters. */
127 : /* -------------------------------------------------------------------- */
128 1 : nXSize = GDALGetRasterBandXSize( hRed );
129 1 : nYSize = GDALGetRasterBandYSize( hRed );
130 :
131 1 : 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 1 : 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 1 : 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 1 : 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 1 : 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 1 : CPLCalloc(GMC_B_LEN * GMC_B_LEN * GMC_B_LEN,sizeof(int));
181 1 : usedboxes = NULL;
182 1 : box_list = freeboxes = (Colorbox *)CPLMalloc(nColors*sizeof (Colorbox));
183 1 : freeboxes[0].next = &freeboxes[1];
184 1 : freeboxes[0].prev = NULL;
185 7 : for (i = 1; i < nColors-1; ++i) {
186 6 : freeboxes[i].next = &freeboxes[i+1];
187 6 : freeboxes[i].prev = &freeboxes[i-1];
188 : }
189 1 : freeboxes[nColors-1].next = NULL;
190 1 : 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 1 : ptr = freeboxes;
202 1 : freeboxes = ptr->next;
203 1 : if (freeboxes)
204 1 : freeboxes->prev = NULL;
205 1 : ptr->next = usedboxes;
206 1 : usedboxes = ptr;
207 1 : if (ptr->next)
208 0 : ptr->next->prev = ptr;
209 :
210 1 : ptr->rmin = ptr->gmin = ptr->bmin = 999;
211 1 : ptr->rmax = ptr->gmax = ptr->bmax = -1;
212 1 : ptr->total = nXSize * nYSize;
213 :
214 : /* -------------------------------------------------------------------- */
215 : /* Collect histogram. */
216 : /* -------------------------------------------------------------------- */
217 1 : pabyRedLine = (GByte *) VSIMalloc(nXSize);
218 1 : pabyGreenLine = (GByte *) VSIMalloc(nXSize);
219 1 : pabyBlueLine = (GByte *) VSIMalloc(nXSize);
220 :
221 1 : 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 51 : for( iLine = 0; iLine < nYSize; iLine++ )
232 : {
233 50 : 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 50 : pabyRedLine, nXSize, 1, GDT_Byte, 0, 0 );
243 : GDALRasterIO( hGreen, GF_Read, 0, iLine, nXSize, 1,
244 50 : pabyGreenLine, nXSize, 1, GDT_Byte, 0, 0 );
245 : GDALRasterIO( hBlue, GF_Read, 0, iLine, nXSize, 1,
246 50 : pabyBlueLine, nXSize, 1, GDT_Byte, 0, 0 );
247 :
248 2550 : for( iPixel = 0; iPixel < nXSize; iPixel++ )
249 : {
250 : int nRed, nGreen, nBlue;
251 :
252 2500 : nRed = pabyRedLine[iPixel] >> COLOR_SHIFT;
253 2500 : nGreen = pabyGreenLine[iPixel] >> COLOR_SHIFT;
254 2500 : nBlue = pabyBlueLine[iPixel] >> COLOR_SHIFT;
255 :
256 2500 : ptr->rmin = MIN(ptr->rmin, nRed);
257 2500 : ptr->gmin = MIN(ptr->gmin, nGreen);
258 2500 : ptr->bmin = MIN(ptr->bmin, nBlue);
259 2500 : ptr->rmax = MAX(ptr->rmax, nRed);
260 2500 : ptr->gmax = MAX(ptr->gmax, nGreen);
261 2500 : ptr->bmax = MAX(ptr->bmax, nBlue);
262 :
263 2500 : histogram[nRed][nGreen][nBlue]++;
264 : }
265 : }
266 :
267 1 : 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 9 : while (freeboxes != NULL) {
279 7 : ptr = largest_box(usedboxes);
280 7 : if (ptr != NULL)
281 7 : splitbox(ptr, histogram, &freeboxes, &usedboxes);
282 : else
283 0 : freeboxes = NULL;
284 : }
285 :
286 : /* ==================================================================== */
287 : /* STEP 4: assign colors to all boxes */
288 : /* ==================================================================== */
289 9 : for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next)
290 : {
291 : GDALColorEntry sEntry;
292 :
293 8 : sEntry.c1 = (GByte) (((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2);
294 8 : sEntry.c2 = (GByte) (((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2);
295 8 : sEntry.c3 = (GByte) (((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2);
296 8 : sEntry.c4 = 255;
297 8 : GDALSetColorEntry( hColorTable, i, &sEntry );
298 : }
299 :
300 : end_and_cleanup:
301 1 : CPLFree( pabyRedLine );
302 1 : CPLFree( pabyGreenLine );
303 1 : CPLFree( pabyBlueLine );
304 :
305 : /* We're done with the boxes now */
306 1 : CPLFree(box_list);
307 1 : freeboxes = usedboxes = NULL;
308 :
309 1 : CPLFree( histogram );
310 :
311 1 : return err;
312 : }
313 :
314 : /************************************************************************/
315 : /* largest_box() */
316 : /************************************************************************/
317 :
318 : static Colorbox *
319 7 : largest_box(Colorbox *usedboxes)
320 : {
321 : Colorbox *p, *b;
322 : int size;
323 :
324 7 : b = NULL;
325 7 : size = -1;
326 35 : for (p = usedboxes; p != NULL; p = p->next)
327 28 : if ((p->rmax > p->rmin || p->gmax > p->gmin ||
328 : p->bmax > p->bmin) && p->total > size)
329 15 : size = (b = p)->total;
330 7 : return (b);
331 : }
332 :
333 : /************************************************************************/
334 : /* splitbox() */
335 : /************************************************************************/
336 : static void
337 7 : splitbox(Colorbox* ptr, int (*histogram)[GMC_B_LEN][GMC_B_LEN],
338 : Colorbox **pfreeboxes, Colorbox **pusedboxes)
339 : {
340 : int hist2[GMC_B_LEN];
341 7 : 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 7 : i = ptr->rmax - ptr->rmin;
355 9 : if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
356 2 : axis = RED;
357 5 : else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
358 4 : axis = GREEN;
359 : else
360 1 : axis = BLUE;
361 : /* get histogram along longest axis */
362 7 : switch (axis) {
363 : case RED:
364 2 : histp = &hist2[ptr->rmin];
365 48 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
366 46 : *histp = 0;
367 1136 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
368 1090 : iptr = &histogram[ir][ig][ptr->bmin];
369 24630 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
370 23540 : *histp += *iptr++;
371 : }
372 46 : histp++;
373 : }
374 2 : first = ptr->rmin;
375 2 : last = ptr->rmax;
376 2 : break;
377 : case GREEN:
378 4 : histp = &hist2[ptr->gmin];
379 58 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
380 54 : *histp = 0;
381 489 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
382 435 : iptr = &histogram[ir][ig][ptr->bmin];
383 3825 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
384 3390 : *histp += *iptr++;
385 : }
386 54 : histp++;
387 : }
388 4 : first = ptr->gmin;
389 4 : last = ptr->gmax;
390 4 : break;
391 : case BLUE:
392 1 : histp = &hist2[ptr->bmin];
393 23 : for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) {
394 22 : *histp = 0;
395 418 : for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
396 396 : iptr = &histogram[ir][ptr->gmin][ib];
397 7524 : for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
398 7128 : *histp += *iptr;
399 7128 : iptr += GMC_B_LEN;
400 : }
401 : }
402 22 : histp++;
403 : }
404 1 : first = ptr->bmin;
405 1 : last = ptr->bmax;
406 : break;
407 : }
408 : /* find median point */
409 7 : sum2 = ptr->total / 2;
410 7 : histp = &hist2[first];
411 7 : sum = 0;
412 7 : for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
413 : ;
414 7 : if (i == first)
415 1 : i++;
416 :
417 : /* Create new box, re-allocate points */
418 7 : new_cb = *pfreeboxes;
419 7 : *pfreeboxes = new_cb->next;
420 7 : if (*pfreeboxes)
421 6 : (*pfreeboxes)->prev = NULL;
422 7 : if (*pusedboxes)
423 7 : (*pusedboxes)->prev = new_cb;
424 7 : new_cb->next = *pusedboxes;
425 7 : *pusedboxes = new_cb;
426 :
427 7 : histp = &hist2[first];
428 43 : for (sum1 = 0, j = first; j < i; j++)
429 36 : sum1 += *histp++;
430 93 : for (sum2 = 0, j = i; j <= last; j++)
431 86 : sum2 += *histp++;
432 7 : new_cb->total = sum1;
433 7 : ptr->total = sum2;
434 :
435 7 : new_cb->rmin = ptr->rmin;
436 7 : new_cb->rmax = ptr->rmax;
437 7 : new_cb->gmin = ptr->gmin;
438 7 : new_cb->gmax = ptr->gmax;
439 7 : new_cb->bmin = ptr->bmin;
440 7 : new_cb->bmax = ptr->bmax;
441 7 : switch (axis) {
442 : case RED:
443 2 : new_cb->rmax = i-1;
444 2 : ptr->rmin = i;
445 2 : break;
446 : case GREEN:
447 4 : new_cb->gmax = i-1;
448 4 : ptr->gmin = i;
449 4 : break;
450 : case BLUE:
451 1 : new_cb->bmax = i-1;
452 1 : ptr->bmin = i;
453 : break;
454 : }
455 7 : shrinkbox(new_cb, histogram);
456 7 : shrinkbox(ptr, histogram);
457 7 : }
458 :
459 : /************************************************************************/
460 : /* shrinkbox() */
461 : /************************************************************************/
462 : static void
463 14 : shrinkbox(Colorbox* box, int (*histogram)[GMC_B_LEN][GMC_B_LEN])
464 : {
465 : int *histp, ir, ig, ib;
466 :
467 14 : if (box->rmax > box->rmin) {
468 20 : for (ir = box->rmin; ir <= box->rmax; ++ir)
469 64 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
470 58 : histp = &histogram[ir][ig][box->bmin];
471 643 : for (ib = box->bmin; ib <= box->bmax; ++ib)
472 599 : if (*histp++ != 0) {
473 14 : box->rmin = ir;
474 14 : goto have_rmin;
475 : }
476 : }
477 : have_rmin:
478 14 : if (box->rmax > box->rmin)
479 39 : for (ir = box->rmax; ir >= box->rmin; --ir)
480 403 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
481 378 : histp = &histogram[ir][ig][box->bmin];
482 378 : ib = box->bmin;
483 3215 : for (; ib <= box->bmax; ++ib)
484 2851 : if (*histp++ != 0) {
485 14 : box->rmax = ir;
486 14 : goto have_rmax;
487 : }
488 : }
489 : }
490 : have_rmax:
491 14 : if (box->gmax > box->gmin) {
492 26 : for (ig = box->gmin; ig <= box->gmax; ++ig)
493 266 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
494 252 : histp = &histogram[ir][ig][box->bmin];
495 5485 : for (ib = box->bmin; ib <= box->bmax; ++ib)
496 5245 : if (*histp++ != 0) {
497 12 : box->gmin = ig;
498 12 : goto have_gmin;
499 : }
500 : }
501 : have_gmin:
502 12 : if (box->gmax > box->gmin)
503 33 : for (ig = box->gmax; ig >= box->gmin; --ig)
504 258 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
505 237 : histp = &histogram[ir][ig][box->bmin];
506 237 : ib = box->bmin;
507 4361 : for (; ib <= box->bmax; ++ib)
508 4136 : if (*histp++ != 0) {
509 12 : box->gmax = ig;
510 12 : goto have_gmax;
511 : }
512 : }
513 : }
514 : have_gmax:
515 14 : if (box->bmax > box->bmin) {
516 16 : for (ib = box->bmin; ib <= box->bmax; ++ib)
517 39 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
518 37 : histp = &histogram[ir][box->gmin][ib];
519 407 : for (ig = box->gmin; ig <= box->gmax; ++ig) {
520 384 : if (*histp != 0) {
521 14 : box->bmin = ib;
522 14 : goto have_bmin;
523 : }
524 370 : histp += GMC_B_LEN;
525 : }
526 : }
527 : have_bmin:
528 14 : if (box->bmax > box->bmin)
529 49 : for (ib = box->bmax; ib >= box->bmin; --ib)
530 313 : for (ir = box->rmin; ir <= box->rmax; ++ir) {
531 278 : histp = &histogram[ir][box->gmin][ib];
532 278 : ig = box->gmin;
533 4352 : for (; ig <= box->gmax; ++ig) {
534 4088 : if (*histp != 0) {
535 14 : box->bmax = ib;
536 14 : goto have_bmax;
537 : }
538 4074 : histp += GMC_B_LEN;
539 : }
540 : }
541 : }
542 : have_bmax:
543 : ;
544 14 : }
|
