00001 /****************************************************************************** 00002 * $Id: gdalmediancut_cpp-source.html,v 1.2 2002/12/21 19:13:13 warmerda Exp $ 00003 * 00004 * Project: CIETMap Phase 2 00005 * Purpose: Use median cut algorithm to generate an near-optimal PCT for a 00006 * given RGB image. Implemented as function GDALComputeMedianCutPCT. 00007 * Author: Frank Warmerdam, warmerdam@pobox.com 00008 * 00009 ****************************************************************************** 00010 * Copyright (c) 2001, Frank Warmerdam 00011 * 00012 * Permission is hereby granted, free of charge, to any person obtaining a 00013 * copy of this software and associated documentation files (the "Software"), 00014 * to deal in the Software without restriction, including without limitation 00015 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 00016 * and/or sell copies of the Software, and to permit persons to whom the 00017 * Software is furnished to do so, subject to the following conditions: 00018 * 00019 * The above copyright notice and this permission notice shall be included 00020 * in all copies or substantial portions of the Software. 00021 * 00022 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 00023 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 00024 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 00025 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 00026 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 00027 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 00028 * DEALINGS IN THE SOFTWARE. 00029 ****************************************************************************** 00030 * 00031 * This code was based on the tiffmedian.c code from libtiff (www.libtiff.org) 00032 * which was based on a paper by Paul Heckbert: 00033 * 00034 * "Color Image Quantization for Frame Buffer Display", Paul 00035 * Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307. 00036 * 00037 * $Log: gdalmediancut_cpp-source.html,v $ 00037 * Revision 1.2 2002/12/21 19:13:13 warmerda 00037 * updated 00037 * 00038 * Revision 1.3 2002/04/16 17:48:36 warmerda 00039 * Ensure everything is initialized. 00040 * 00041 * Revision 1.2 2001/07/18 04:43:13 warmerda 00042 * added CPL_CVSID 00043 * 00044 * Revision 1.1 2001/01/22 22:30:59 warmerda 00045 * New 00046 * 00047 */ 00048 00049 #include "gdal_priv.h" 00050 #include "gdal_alg.h" 00051 00052 CPL_CVSID("$Id: gdalmediancut_cpp-source.html,v 1.2 2002/12/21 19:13:13 warmerda Exp $"); 00053 00054 #define MAX_CMAP_SIZE 256 00055 00056 #define COLOR_DEPTH 8 00057 #define MAX_COLOR 256 00058 00059 #define GMC_B_DEPTH 5 /* # bits/pixel to use */ 00060 #define GMC_B_LEN (1L<<GMC_B_DEPTH) 00061 00062 #define C_DEPTH 2 00063 #define C_LEN (1L<<C_DEPTH) /* # cells/color to use */ 00064 00065 #define COLOR_SHIFT (COLOR_DEPTH-GMC_B_DEPTH) 00066 00067 typedef struct colorbox { 00068 struct colorbox *next, *prev; 00069 int rmin, rmax; 00070 int gmin, gmax; 00071 int bmin, bmax; 00072 int total; 00073 } Colorbox; 00074 00075 static int num_colors; 00076 static int (*histogram)[GMC_B_LEN][GMC_B_LEN]; 00077 static Colorbox *freeboxes; 00078 static Colorbox *usedboxes; 00079 00080 static void splitbox(Colorbox*); 00081 static void shrinkbox(Colorbox*); 00082 static Colorbox* largest_box(void); 00083 00084 /************************************************************************/ 00085 /* GDALComputeMedianCutPCT() */ 00086 /************************************************************************/ 00087 00088 int GDALComputeMedianCutPCT( GDALRasterBandH hRed, 00089 GDALRasterBandH hGreen, 00090 GDALRasterBandH hBlue, 00091 int (*pfnIncludePixel)(int,int,void*), 00092 int nColors, 00093 GDALColorTableH hColorTable, 00094 GDALProgressFunc pfnProgress, 00095 void * pProgressArg ) 00096 00097 { 00098 int nXSize, nYSize; 00099 00100 /* -------------------------------------------------------------------- */ 00101 /* Validate parameters. */ 00102 /* -------------------------------------------------------------------- */ 00103 nXSize = GDALGetRasterBandXSize( hRed ); 00104 nYSize = GDALGetRasterBandYSize( hRed ); 00105 00106 if( GDALGetRasterBandXSize( hGreen ) != nXSize 00107 || GDALGetRasterBandYSize( hGreen ) != nYSize 00108 || GDALGetRasterBandXSize( hBlue ) != nXSize 00109 || GDALGetRasterBandYSize( hBlue ) != nYSize ) 00110 { 00111 CPLError( CE_Failure, CPLE_IllegalArg, 00112 "Green or blue band doesn't match size of red band.\n" ); 00113 00114 return CE_Failure; 00115 } 00116 00117 if( pfnIncludePixel != NULL ) 00118 { 00119 CPLError( CE_Failure, CPLE_IllegalArg, 00120 "GDALComputeMedianCutPCT() doesn't currently support " 00121 " pfnIncludePixel function." ); 00122 00123 return CE_Failure; 00124 } 00125 00126 if( pfnProgress == NULL ) 00127 pfnProgress = GDALDummyProgress; 00128 00129 histogram = (int (*)[GMC_B_LEN][GMC_B_LEN]) 00130 CPLCalloc(GMC_B_LEN * GMC_B_LEN * GMC_B_LEN,sizeof(int)); 00131 00132 /* ==================================================================== */ 00133 /* STEP 1: crate empty boxes. */ 00134 /* ==================================================================== */ 00135 int i; 00136 Colorbox *box_list, *ptr; 00137 00138 num_colors = nColors; 00139 usedboxes = NULL; 00140 box_list = freeboxes = (Colorbox *)CPLMalloc(num_colors*sizeof (Colorbox)); 00141 freeboxes[0].next = &freeboxes[1]; 00142 freeboxes[0].prev = NULL; 00143 for (i = 1; i < num_colors-1; ++i) { 00144 freeboxes[i].next = &freeboxes[i+1]; 00145 freeboxes[i].prev = &freeboxes[i-1]; 00146 } 00147 freeboxes[num_colors-1].next = NULL; 00148 freeboxes[num_colors-1].prev = &freeboxes[num_colors-2]; 00149 00150 /* ==================================================================== */ 00151 /* Build histogram. */ 00152 /* ==================================================================== */ 00153 GByte *pabyRedLine, *pabyGreenLine, *pabyBlueLine; 00154 int iLine, iPixel; 00155 00156 /* -------------------------------------------------------------------- */ 00157 /* Initialize the box datastructures. */ 00158 /* -------------------------------------------------------------------- */ 00159 ptr = freeboxes; 00160 freeboxes = ptr->next; 00161 if (freeboxes) 00162 freeboxes->prev = NULL; 00163 ptr->next = usedboxes; 00164 usedboxes = ptr; 00165 if (ptr->next) 00166 ptr->next->prev = ptr; 00167 00168 ptr->rmin = ptr->gmin = ptr->bmin = 999; 00169 ptr->rmax = ptr->gmax = ptr->bmax = -1; 00170 ptr->total = nXSize * nYSize; 00171 00172 memset( histogram, 0, sizeof(int) * GMC_B_LEN * GMC_B_LEN * GMC_B_LEN ); 00173 00174 /* -------------------------------------------------------------------- */ 00175 /* Collect histogram. */ 00176 /* -------------------------------------------------------------------- */ 00177 pabyRedLine = (GByte *) CPLMalloc(nXSize); 00178 pabyGreenLine = (GByte *) CPLMalloc(nXSize); 00179 pabyBlueLine = (GByte *) CPLMalloc(nXSize); 00180 00181 for( iLine = 0; iLine < nYSize; iLine++ ) 00182 { 00183 if( !pfnProgress( iLine / (double) nYSize, 00184 "Generating Histogram", pProgressArg ) ) 00185 { 00186 CPLFree( pabyRedLine ); 00187 CPLFree( pabyGreenLine ); 00188 CPLFree( pabyBlueLine ); 00189 00190 CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" ); 00191 return CE_Failure; 00192 } 00193 00194 GDALRasterIO( hRed, GF_Read, 0, iLine, nXSize, 1, 00195 pabyRedLine, nXSize, 1, GDT_Byte, 0, 0 ); 00196 GDALRasterIO( hGreen, GF_Read, 0, iLine, nXSize, 1, 00197 pabyGreenLine, nXSize, 1, GDT_Byte, 0, 0 ); 00198 GDALRasterIO( hBlue, GF_Read, 0, iLine, nXSize, 1, 00199 pabyBlueLine, nXSize, 1, GDT_Byte, 0, 0 ); 00200 00201 for( iPixel = 0; iPixel < nXSize; iPixel++ ) 00202 { 00203 int nRed, nGreen, nBlue; 00204 00205 nRed = pabyRedLine[iPixel] >> COLOR_SHIFT; 00206 nGreen = pabyGreenLine[iPixel] >> COLOR_SHIFT; 00207 nBlue = pabyBlueLine[iPixel] >> COLOR_SHIFT; 00208 00209 ptr->rmin = MIN(ptr->rmin, nRed); 00210 ptr->gmin = MIN(ptr->gmin, nGreen); 00211 ptr->bmin = MIN(ptr->bmin, nBlue); 00212 ptr->rmax = MAX(ptr->rmax, nRed); 00213 ptr->gmax = MAX(ptr->gmax, nGreen); 00214 ptr->bmax = MAX(ptr->bmax, nBlue); 00215 00216 histogram[nRed][nGreen][nBlue]++; 00217 } 00218 } 00219 00220 CPLFree( pabyRedLine ); 00221 CPLFree( pabyGreenLine ); 00222 CPLFree( pabyBlueLine ); 00223 00224 if( !pfnProgress( 1.0, "Generating Histogram", pProgressArg ) ) 00225 { 00226 CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" ); 00227 return CE_Failure; 00228 } 00229 00230 /* ==================================================================== */ 00231 /* STEP 3: continually subdivide boxes until no more free */ 00232 /* boxes remain or until all colors assigned. */ 00233 /* ==================================================================== */ 00234 while (freeboxes != NULL) { 00235 ptr = largest_box(); 00236 if (ptr != NULL) 00237 splitbox(ptr); 00238 else 00239 freeboxes = NULL; 00240 } 00241 00242 /* ==================================================================== */ 00243 /* STEP 4: assign colors to all boxes */ 00244 /* ==================================================================== */ 00245 for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next) 00246 { 00247 GDALColorEntry sEntry; 00248 00249 sEntry.c1 = ((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2; 00250 sEntry.c2 = ((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2; 00251 sEntry.c3 = ((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2; 00252 GDALSetColorEntry( hColorTable, i, &sEntry ); 00253 } 00254 00255 /* We're done with the boxes now */ 00256 CPLFree(box_list); 00257 freeboxes = usedboxes = NULL; 00258 00259 CPLFree( histogram ); 00260 00261 return CE_None; 00262 } 00263 00264 /************************************************************************/ 00265 /* largest_box() */ 00266 /************************************************************************/ 00267 00268 static Colorbox * 00269 largest_box(void) 00270 { 00271 register Colorbox *p, *b; 00272 register int size; 00273 00274 b = NULL; 00275 size = -1; 00276 for (p = usedboxes; p != NULL; p = p->next) 00277 if ((p->rmax > p->rmin || p->gmax > p->gmin || 00278 p->bmax > p->bmin) && p->total > size) 00279 size = (b = p)->total; 00280 return (b); 00281 } 00282 00283 /************************************************************************/ 00284 /* splitbox() */ 00285 /************************************************************************/ 00286 static void 00287 splitbox(Colorbox* ptr) 00288 { 00289 int hist2[GMC_B_LEN]; 00290 int first=0, last=0; 00291 register Colorbox *new_cb; 00292 register int *iptr, *histp; 00293 register int i, j; 00294 register int ir,ig,ib; 00295 register int sum, sum1, sum2; 00296 enum { RED, GREEN, BLUE } axis; 00297 00298 /* 00299 * See which axis is the largest, do a histogram along that 00300 * axis. Split at median point. Contract both new boxes to 00301 * fit points and return 00302 */ 00303 i = ptr->rmax - ptr->rmin; 00304 if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin) 00305 axis = RED; 00306 else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin) 00307 axis = GREEN; 00308 else 00309 axis = BLUE; 00310 /* get histogram along longest axis */ 00311 switch (axis) { 00312 case RED: 00313 histp = &hist2[ptr->rmin]; 00314 for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) { 00315 *histp = 0; 00316 for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) { 00317 iptr = &histogram[ir][ig][ptr->bmin]; 00318 for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) 00319 *histp += *iptr++; 00320 } 00321 histp++; 00322 } 00323 first = ptr->rmin; 00324 last = ptr->rmax; 00325 break; 00326 case GREEN: 00327 histp = &hist2[ptr->gmin]; 00328 for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) { 00329 *histp = 0; 00330 for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) { 00331 iptr = &histogram[ir][ig][ptr->bmin]; 00332 for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) 00333 *histp += *iptr++; 00334 } 00335 histp++; 00336 } 00337 first = ptr->gmin; 00338 last = ptr->gmax; 00339 break; 00340 case BLUE: 00341 histp = &hist2[ptr->bmin]; 00342 for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) { 00343 *histp = 0; 00344 for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) { 00345 iptr = &histogram[ir][ptr->gmin][ib]; 00346 for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) { 00347 *histp += *iptr; 00348 iptr += GMC_B_LEN; 00349 } 00350 } 00351 histp++; 00352 } 00353 first = ptr->bmin; 00354 last = ptr->bmax; 00355 break; 00356 } 00357 /* find median point */ 00358 sum2 = ptr->total / 2; 00359 histp = &hist2[first]; 00360 sum = 0; 00361 for (i = first; i <= last && (sum += *histp++) < sum2; ++i) 00362 ; 00363 if (i == first) 00364 i++; 00365 00366 /* Create new box, re-allocate points */ 00367 new_cb = freeboxes; 00368 freeboxes = new_cb->next; 00369 if (freeboxes) 00370 freeboxes->prev = NULL; 00371 if (usedboxes) 00372 usedboxes->prev = new_cb; 00373 new_cb->next = usedboxes; 00374 usedboxes = new_cb; 00375 00376 histp = &hist2[first]; 00377 for (sum1 = 0, j = first; j < i; j++) 00378 sum1 += *histp++; 00379 for (sum2 = 0, j = i; j <= last; j++) 00380 sum2 += *histp++; 00381 new_cb->total = sum1; 00382 ptr->total = sum2; 00383 00384 new_cb->rmin = ptr->rmin; 00385 new_cb->rmax = ptr->rmax; 00386 new_cb->gmin = ptr->gmin; 00387 new_cb->gmax = ptr->gmax; 00388 new_cb->bmin = ptr->bmin; 00389 new_cb->bmax = ptr->bmax; 00390 switch (axis) { 00391 case RED: 00392 new_cb->rmax = i-1; 00393 ptr->rmin = i; 00394 break; 00395 case GREEN: 00396 new_cb->gmax = i-1; 00397 ptr->gmin = i; 00398 break; 00399 case BLUE: 00400 new_cb->bmax = i-1; 00401 ptr->bmin = i; 00402 break; 00403 } 00404 shrinkbox(new_cb); 00405 shrinkbox(ptr); 00406 } 00407 00408 /************************************************************************/ 00409 /* shrinkbox() */ 00410 /************************************************************************/ 00411 static void 00412 shrinkbox(Colorbox* box) 00413 { 00414 register int *histp, ir, ig, ib; 00415 00416 if (box->rmax > box->rmin) { 00417 for (ir = box->rmin; ir <= box->rmax; ++ir) 00418 for (ig = box->gmin; ig <= box->gmax; ++ig) { 00419 histp = &histogram[ir][ig][box->bmin]; 00420 for (ib = box->bmin; ib <= box->bmax; ++ib) 00421 if (*histp++ != 0) { 00422 box->rmin = ir; 00423 goto have_rmin; 00424 } 00425 } 00426 have_rmin: 00427 if (box->rmax > box->rmin) 00428 for (ir = box->rmax; ir >= box->rmin; --ir) 00429 for (ig = box->gmin; ig <= box->gmax; ++ig) { 00430 histp = &histogram[ir][ig][box->bmin]; 00431 ib = box->bmin; 00432 for (; ib <= box->bmax; ++ib) 00433 if (*histp++ != 0) { 00434 box->rmax = ir; 00435 goto have_rmax; 00436 } 00437 } 00438 } 00439 have_rmax: 00440 if (box->gmax > box->gmin) { 00441 for (ig = box->gmin; ig <= box->gmax; ++ig) 00442 for (ir = box->rmin; ir <= box->rmax; ++ir) { 00443 histp = &histogram[ir][ig][box->bmin]; 00444 for (ib = box->bmin; ib <= box->bmax; ++ib) 00445 if (*histp++ != 0) { 00446 box->gmin = ig; 00447 goto have_gmin; 00448 } 00449 } 00450 have_gmin: 00451 if (box->gmax > box->gmin) 00452 for (ig = box->gmax; ig >= box->gmin; --ig) 00453 for (ir = box->rmin; ir <= box->rmax; ++ir) { 00454 histp = &histogram[ir][ig][box->bmin]; 00455 ib = box->bmin; 00456 for (; ib <= box->bmax; ++ib) 00457 if (*histp++ != 0) { 00458 box->gmax = ig; 00459 goto have_gmax; 00460 } 00461 } 00462 } 00463 have_gmax: 00464 if (box->bmax > box->bmin) { 00465 for (ib = box->bmin; ib <= box->bmax; ++ib) 00466 for (ir = box->rmin; ir <= box->rmax; ++ir) { 00467 histp = &histogram[ir][box->gmin][ib]; 00468 for (ig = box->gmin; ig <= box->gmax; ++ig) { 00469 if (*histp != 0) { 00470 box->bmin = ib; 00471 goto have_bmin; 00472 } 00473 histp += GMC_B_LEN; 00474 } 00475 } 00476 have_bmin: 00477 if (box->bmax > box->bmin) 00478 for (ib = box->bmax; ib >= box->bmin; --ib) 00479 for (ir = box->rmin; ir <= box->rmax; ++ir) { 00480 histp = &histogram[ir][box->gmin][ib]; 00481 ig = box->gmin; 00482 for (; ig <= box->gmax; ++ig) { 00483 if (*histp != 0) { 00484 box->bmax = ib; 00485 goto have_bmax; 00486 } 00487 histp += GMC_B_LEN; 00488 } 00489 } 00490 } 00491 have_bmax: 00492 ; 00493 }