Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members   Related Pages  

gdalmediancut.cpp

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 }

Generated at Sat Dec 21 14:01:59 2002 for GDAL by doxygen1.2.3-20001105 written by Dimitri van Heesch, © 1997-2000