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 }