LCOV - code coverage report
Current view: top level - alg - gdalmediancut.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 525 562 93.4 %
Date: 2025-02-20 10:14:44 Functions: 13 18 72.2 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  CIETMap Phase 2
       4             :  * Purpose:  Use median cut algorithm to generate an near-optimal PCT for a
       5             :  *           given RGB image.  Implemented as function GDALComputeMedianCutPCT.
       6             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2001, Frank Warmerdam
      10             :  * Copyright (c) 2007-2010, Even Rouault <even dot rouault at spatialys.com>
      11             :  *
      12             :  * SPDX-License-Identifier: MIT
      13             :  ******************************************************************************
      14             :  *
      15             :  * This code was based on the tiffmedian.c code from libtiff (www.libtiff.org)
      16             :  * which was based on a paper by Paul Heckbert:
      17             :  *
      18             :  *      "Color  Image Quantization for Frame Buffer Display", Paul
      19             :  *      Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307.
      20             :  *
      21             :  */
      22             : 
      23             : #include "cpl_port.h"
      24             : #include "gdal_alg.h"
      25             : #include "gdal_alg_priv.h"
      26             : 
      27             : #include <climits>
      28             : #include <cstring>
      29             : 
      30             : #include <algorithm>
      31             : #include <limits>
      32             : 
      33             : #include "cpl_conv.h"
      34             : #include "cpl_error.h"
      35             : #include "cpl_float.h"
      36             : #include "cpl_progress.h"
      37             : #include "cpl_vsi.h"
      38             : #include "gdal.h"
      39             : #include "gdal_priv.h"
      40             : 
      41      287563 : template <typename T> static T *HISTOGRAM(T *h, int n, int r, int g, int b)
      42             : {
      43      287563 :     const int index = (r * n + g) * n + b;
      44      287563 :     return &h[index];
      45             : }
      46             : 
      47             : #ifndef MAKE_COLOR_CODE_defined
      48             : #define MAKE_COLOR_CODE_defined
      49             : 
      50    17939300 : static int MAKE_COLOR_CODE(int r, int g, int b)
      51             : {
      52    17939300 :     return r | (g << 8) | (b << 16);
      53             : }
      54             : #endif
      55             : 
      56             : // NOTE: If changing the size of this structure, edit
      57             : // MEDIAN_CUT_AND_DITHER_BUFFER_SIZE_65536 in gdal_alg_priv.h and take into
      58             : // account ColorIndex in gdaldither.cpp.
      59             : typedef struct
      60             : {
      61             :     GUInt32 nColorCode;
      62             :     int nCount;
      63             :     GUInt32 nColorCode2;
      64             :     int nCount2;
      65             :     GUInt32 nColorCode3;
      66             :     int nCount3;
      67             : } HashHistogram;
      68             : 
      69             : typedef struct colorbox
      70             : {
      71             :     struct colorbox *next, *prev;
      72             :     int rmin, rmax;
      73             :     int gmin, gmax;
      74             :     int bmin, bmax;
      75             :     GUIntBig total;
      76             : } Colorbox;
      77             : 
      78             : template <class T>
      79             : static void splitbox(Colorbox *ptr, const T *histogram,
      80             :                      const HashHistogram *psHashHistogram, int nCLevels,
      81             :                      Colorbox **pfreeboxes, Colorbox **pusedboxes,
      82             :                      GByte *pabyRedBand, GByte *pabyGreenBand,
      83             :                      GByte *pabyBlueBand, T nPixels);
      84             : 
      85             : template <class T>
      86             : static void shrinkbox(Colorbox *box, const T *histogram, int nCLevels);
      87             : 
      88             : static Colorbox *largest_box(Colorbox *usedboxes);
      89             : 
      90             : /************************************************************************/
      91             : /*                      GDALComputeMedianCutPCT()                       */
      92             : /************************************************************************/
      93             : 
      94             : /**
      95             :  * Compute optimal PCT for RGB image.
      96             :  *
      97             :  * This function implements a median cut algorithm to compute an "optimal"
      98             :  * pseudocolor table for representing an input RGB image.  This PCT could
      99             :  * then be used with GDALDitherRGB2PCT() to convert a 24bit RGB image into
     100             :  * an eightbit pseudo-colored image.
     101             :  *
     102             :  * This code was based on the tiffmedian.c code from libtiff (www.libtiff.org)
     103             :  * which was based on a paper by Paul Heckbert:
     104             :  *
     105             :  * \verbatim
     106             :  *   "Color  Image Quantization for Frame Buffer Display", Paul
     107             :  *   Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307.
     108             :  * \endverbatim
     109             :  *
     110             :  * The red, green and blue input bands do not necessarily need to come
     111             :  * from the same file, but they must be the same width and height.  They will
     112             :  * be clipped to 8bit during reading, so non-eight bit bands are generally
     113             :  * inappropriate.
     114             :  *
     115             :  * @param hRed Red input band.
     116             :  * @param hGreen Green input band.
     117             :  * @param hBlue Blue input band.
     118             :  * @param pfnIncludePixel function used to test which pixels should be included
     119             :  * in the analysis.  At this time this argument is ignored and all pixels are
     120             :  * utilized.  This should normally be NULL.
     121             :  * @param nColors the desired number of colors to be returned (2-256).
     122             :  * @param hColorTable the colors will be returned in this color table object.
     123             :  * @param pfnProgress callback for reporting algorithm progress matching the
     124             :  * GDALProgressFunc() semantics.  May be NULL.
     125             :  * @param pProgressArg callback argument passed to pfnProgress.
     126             :  *
     127             :  * @return returns CE_None on success or CE_Failure if an error occurs.
     128             :  */
     129             : 
     130           3 : extern "C" int CPL_STDCALL GDALComputeMedianCutPCT(
     131             :     GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
     132             :     int (*pfnIncludePixel)(int, int, void *), int nColors,
     133             :     GDALColorTableH hColorTable, GDALProgressFunc pfnProgress,
     134             :     void *pProgressArg)
     135             : 
     136             : {
     137           3 :     VALIDATE_POINTER1(hRed, "GDALComputeMedianCutPCT", CE_Failure);
     138           3 :     const int nXSize = GDALGetRasterBandXSize(hRed);
     139           3 :     const int nYSize = GDALGetRasterBandYSize(hRed);
     140           3 :     if (nYSize == 0)
     141           0 :         return CE_Failure;
     142           3 :     if (static_cast<GUInt32>(nXSize) <
     143           3 :         std::numeric_limits<GUInt32>::max() / static_cast<GUInt32>(nYSize))
     144             :     {
     145           3 :         return GDALComputeMedianCutPCTInternal(
     146             :             hRed, hGreen, hBlue, nullptr, nullptr, nullptr, pfnIncludePixel,
     147             :             nColors, 5, static_cast<GUInt32 *>(nullptr), hColorTable,
     148           3 :             pfnProgress, pProgressArg);
     149             :     }
     150             :     else
     151             :     {
     152           0 :         return GDALComputeMedianCutPCTInternal(
     153             :             hRed, hGreen, hBlue, nullptr, nullptr, nullptr, pfnIncludePixel,
     154             :             nColors, 5, static_cast<GUIntBig *>(nullptr), hColorTable,
     155           0 :             pfnProgress, pProgressArg);
     156             :     }
     157             : }
     158             : 
     159             : #ifndef IsColorCodeSet_defined
     160             : #define IsColorCodeSet_defined
     161             : 
     162    18835800 : static inline bool IsColorCodeSet(GUInt32 nColorCode)
     163             : {
     164    18835800 :     return (nColorCode >> 31) == 0;
     165             : }
     166             : #endif
     167             : 
     168    17549600 : static inline int FindColorCount(const HashHistogram *psHashHistogram,
     169             :                                  GUInt32 nColorCode)
     170             : {
     171             : 
     172    17549600 :     GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
     173             :     while (true)
     174             :     {
     175    17550200 :         if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode))
     176             :         {
     177    16219100 :             return 0;
     178             :         }
     179     1331150 :         if (psHashHistogram[nIdx].nColorCode == nColorCode)
     180             :         {
     181      141464 :             return psHashHistogram[nIdx].nCount;
     182             :         }
     183     1189680 :         if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode2))
     184             :         {
     185     1151920 :             return 0;
     186             :         }
     187       37761 :         if (psHashHistogram[nIdx].nColorCode2 == nColorCode)
     188             :         {
     189        4030 :             return psHashHistogram[nIdx].nCount2;
     190             :         }
     191       33731 :         if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode3))
     192             :         {
     193       33078 :             return 0;
     194             :         }
     195         653 :         if (psHashHistogram[nIdx].nColorCode3 == nColorCode)
     196             :         {
     197          59 :             return psHashHistogram[nIdx].nCount3;
     198             :         }
     199             : 
     200           0 :         do
     201             :         {
     202         594 :             nIdx += 257;
     203         594 :             if (nIdx >= PRIME_FOR_65536)
     204           0 :                 nIdx -= PRIME_FOR_65536;
     205         594 :         } while (IsColorCodeSet(psHashHistogram[nIdx].nColorCode) &&
     206         260 :                  psHashHistogram[nIdx].nColorCode != nColorCode &&
     207         130 :                  IsColorCodeSet(psHashHistogram[nIdx].nColorCode2) &&
     208           0 :                  psHashHistogram[nIdx].nColorCode2 != nColorCode &&
     209         724 :                  IsColorCodeSet(psHashHistogram[nIdx].nColorCode3) &&
     210           0 :                  psHashHistogram[nIdx].nColorCode3 != nColorCode);
     211             :     }
     212             : }
     213             : 
     214      389644 : static inline int *FindAndInsertColorCount(HashHistogram *psHashHistogram,
     215             :                                            GUInt32 nColorCode)
     216             : {
     217      389644 :     GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
     218             :     while (true)
     219             :     {
     220      389644 :         if (psHashHistogram[nIdx].nColorCode == nColorCode)
     221             :         {
     222      329957 :             return &(psHashHistogram[nIdx].nCount);
     223             :         }
     224       59687 :         if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode))
     225             :         {
     226       56306 :             psHashHistogram[nIdx].nColorCode = nColorCode;
     227       56306 :             psHashHistogram[nIdx].nCount = 0;
     228       56306 :             return &(psHashHistogram[nIdx].nCount);
     229             :         }
     230        3381 :         if (psHashHistogram[nIdx].nColorCode2 == nColorCode)
     231             :         {
     232        1661 :             return &(psHashHistogram[nIdx].nCount2);
     233             :         }
     234        1720 :         if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode2))
     235             :         {
     236        1687 :             psHashHistogram[nIdx].nColorCode2 = nColorCode;
     237        1687 :             psHashHistogram[nIdx].nCount2 = 0;
     238        1687 :             return &(psHashHistogram[nIdx].nCount2);
     239             :         }
     240          33 :         if (psHashHistogram[nIdx].nColorCode3 == nColorCode)
     241             :         {
     242           5 :             return &(psHashHistogram[nIdx].nCount3);
     243             :         }
     244          28 :         if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode3))
     245             :         {
     246          28 :             psHashHistogram[nIdx].nColorCode3 = nColorCode;
     247          28 :             psHashHistogram[nIdx].nCount3 = 0;
     248          28 :             return &(psHashHistogram[nIdx].nCount3);
     249             :         }
     250             : 
     251           0 :         do
     252             :         {
     253           0 :             nIdx += 257;
     254           0 :             if (nIdx >= PRIME_FOR_65536)
     255           0 :                 nIdx -= PRIME_FOR_65536;
     256           0 :         } while (IsColorCodeSet(psHashHistogram[nIdx].nColorCode) &&
     257           0 :                  psHashHistogram[nIdx].nColorCode != nColorCode &&
     258           0 :                  IsColorCodeSet(psHashHistogram[nIdx].nColorCode2) &&
     259           0 :                  psHashHistogram[nIdx].nColorCode2 != nColorCode &&
     260           0 :                  IsColorCodeSet(psHashHistogram[nIdx].nColorCode3) &&
     261           0 :                  psHashHistogram[nIdx].nColorCode3 != nColorCode);
     262             :     }
     263             : }
     264             : 
     265             : template <class T>
     266          13 : int GDALComputeMedianCutPCTInternal(
     267             :     GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
     268             :     GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
     269             :     int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
     270             :     T *panHistogram,  // NULL, or >= size (1<<nBits)^3 * sizeof(T) bytes.
     271             :     GDALColorTableH hColorTable, GDALProgressFunc pfnProgress,
     272             :     void *pProgressArg)
     273             : 
     274             : {
     275          13 :     VALIDATE_POINTER1(hRed, "GDALComputeMedianCutPCT", CE_Failure);
     276          13 :     VALIDATE_POINTER1(hGreen, "GDALComputeMedianCutPCT", CE_Failure);
     277          13 :     VALIDATE_POINTER1(hBlue, "GDALComputeMedianCutPCT", CE_Failure);
     278             : 
     279          13 :     CPLErr err = CE_None;
     280             : 
     281             :     /* -------------------------------------------------------------------- */
     282             :     /*      Validate parameters.                                            */
     283             :     /* -------------------------------------------------------------------- */
     284          13 :     const int nXSize = GDALGetRasterBandXSize(hRed);
     285          13 :     const int nYSize = GDALGetRasterBandYSize(hRed);
     286             : 
     287          13 :     if (GDALGetRasterBandXSize(hGreen) != nXSize ||
     288          13 :         GDALGetRasterBandYSize(hGreen) != nYSize ||
     289          39 :         GDALGetRasterBandXSize(hBlue) != nXSize ||
     290          13 :         GDALGetRasterBandYSize(hBlue) != nYSize)
     291             :     {
     292           0 :         CPLError(CE_Failure, CPLE_IllegalArg,
     293             :                  "Green or blue band doesn't match size of red band.");
     294             : 
     295           0 :         return CE_Failure;
     296             :     }
     297             : 
     298          13 :     if (pfnIncludePixel != nullptr)
     299             :     {
     300           0 :         CPLError(CE_Failure, CPLE_IllegalArg,
     301             :                  "GDALComputeMedianCutPCT() doesn't currently support "
     302             :                  "pfnIncludePixel function.");
     303             : 
     304           0 :         return CE_Failure;
     305             :     }
     306             : 
     307          13 :     if (nColors <= 0)
     308             :     {
     309           0 :         CPLError(CE_Failure, CPLE_IllegalArg,
     310             :                  "GDALComputeMedianCutPCT(): "
     311             :                  "nColors must be strictly greater than 1.");
     312             : 
     313           0 :         return CE_Failure;
     314             :     }
     315             : 
     316          13 :     if (nColors > 256)
     317             :     {
     318           0 :         CPLError(CE_Failure, CPLE_IllegalArg,
     319             :                  "GDALComputeMedianCutPCT(): "
     320             :                  "nColors must be lesser than or equal to 256.");
     321             : 
     322           0 :         return CE_Failure;
     323             :     }
     324             : 
     325          13 :     if (pfnProgress == nullptr)
     326          11 :         pfnProgress = GDALDummyProgress;
     327             : 
     328             :     /* ==================================================================== */
     329             :     /*      STEP 1: create empty boxes.                                     */
     330             :     /* ==================================================================== */
     331          13 :     if (static_cast<GUInt32>(nXSize) >
     332          13 :         cpl::NumericLimits<T>::max() / static_cast<GUInt32>(nYSize))
     333             :     {
     334           0 :         CPLError(CE_Warning, CPLE_AppDefined,
     335             :                  "GDALComputeMedianCutPCTInternal() not called "
     336             :                  "with large enough type");
     337             :     }
     338             : 
     339          13 :     T nPixels = 0;
     340          10 :     if (nBits == 8 && pabyRedBand != nullptr && pabyGreenBand != nullptr &&
     341          23 :         pabyBlueBand != nullptr &&
     342          10 :         static_cast<GUInt32>(nXSize) <=
     343          10 :             cpl::NumericLimits<T>::max() / static_cast<GUInt32>(nYSize))
     344             :     {
     345          10 :         nPixels = static_cast<T>(nXSize) * static_cast<T>(nYSize);
     346             :     }
     347             : 
     348          13 :     const int nCLevels = 1 << nBits;
     349          13 :     const int nCLevelsCube = nCLevels * nCLevels * nCLevels;
     350          13 :     T *histogram = nullptr;
     351          13 :     HashHistogram *psHashHistogram = nullptr;
     352          13 :     if (panHistogram)
     353             :     {
     354          10 :         if (nBits == 8 && static_cast<GUIntBig>(nXSize) * nYSize <= 65536)
     355             :         {
     356             :             // If the image is small enough, then the number of colors
     357             :             // will be limited and using a hashmap, rather than a full table
     358             :             // will be more efficient.
     359           9 :             histogram = nullptr;
     360           9 :             psHashHistogram = reinterpret_cast<HashHistogram *>(panHistogram);
     361           9 :             memset(psHashHistogram, 0xFF,
     362             :                    sizeof(HashHistogram) * PRIME_FOR_65536);
     363             :         }
     364             :         else
     365             :         {
     366           1 :             histogram = panHistogram;
     367           1 :             memset(histogram, 0, nCLevelsCube * sizeof(T));
     368             :         }
     369             :     }
     370             :     else
     371             :     {
     372             :         histogram =
     373           3 :             static_cast<T *>(VSI_CALLOC_VERBOSE(nCLevelsCube, sizeof(T)));
     374           3 :         if (histogram == nullptr)
     375             :         {
     376           0 :             return CE_Failure;
     377             :         }
     378             :     }
     379             :     Colorbox *box_list =
     380          13 :         static_cast<Colorbox *>(CPLMalloc(nColors * sizeof(Colorbox)));
     381          13 :     Colorbox *freeboxes = box_list;
     382          13 :     freeboxes[0].next = &freeboxes[1];
     383          13 :     freeboxes[0].prev = nullptr;
     384        2827 :     for (int i = 1; i < nColors - 1; ++i)
     385             :     {
     386        2814 :         freeboxes[i].next = &freeboxes[i + 1];
     387        2814 :         freeboxes[i].prev = &freeboxes[i - 1];
     388             :     }
     389          13 :     freeboxes[nColors - 1].next = nullptr;
     390          13 :     freeboxes[nColors - 1].prev = &freeboxes[nColors - 2];
     391             : 
     392             :     /* ==================================================================== */
     393             :     /*      Build histogram.                                                */
     394             :     /* ==================================================================== */
     395             : 
     396             :     /* -------------------------------------------------------------------- */
     397             :     /*      Initialize the box datastructures.                              */
     398             :     /* -------------------------------------------------------------------- */
     399          13 :     Colorbox *freeboxes_before = freeboxes;
     400          13 :     freeboxes = freeboxes_before->next;
     401          13 :     if (freeboxes)
     402          13 :         freeboxes->prev = nullptr;
     403             : 
     404          13 :     Colorbox *usedboxes = freeboxes_before;
     405          13 :     usedboxes->next = nullptr;
     406          13 :     usedboxes->rmin = 999;
     407          13 :     usedboxes->gmin = 999;
     408          13 :     usedboxes->bmin = 999;
     409          13 :     usedboxes->rmax = -1;
     410          13 :     usedboxes->gmax = -1;
     411          13 :     usedboxes->bmax = -1;
     412          13 :     usedboxes->total =
     413          13 :         static_cast<GUIntBig>(nXSize) * static_cast<GUIntBig>(nYSize);
     414             : 
     415             :     /* -------------------------------------------------------------------- */
     416             :     /*      Collect histogram.                                              */
     417             :     /* -------------------------------------------------------------------- */
     418             : 
     419             :     // TODO(schwehr): Move these closer to usage after removing gotos.
     420          13 :     const int nColorShift = 8 - nBits;
     421          13 :     int nColorCounter = 0;
     422          13 :     GByte anRed[256] = {};
     423          13 :     GByte anGreen[256] = {};
     424          13 :     GByte anBlue[256] = {};
     425             : 
     426          13 :     GByte *pabyRedLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
     427          13 :     GByte *pabyGreenLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
     428          13 :     GByte *pabyBlueLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
     429             : 
     430          13 :     if (pabyRedLine == nullptr || pabyGreenLine == nullptr ||
     431             :         pabyBlueLine == nullptr)
     432             :     {
     433           0 :         err = CE_Failure;
     434           0 :         goto end_and_cleanup;
     435             :     }
     436             : 
     437        2137 :     for (int iLine = 0; iLine < nYSize; iLine++)
     438             :     {
     439        2124 :         if (!pfnProgress(iLine / static_cast<double>(nYSize),
     440             :                          "Generating Histogram", pProgressArg))
     441             :         {
     442           0 :             CPLError(CE_Failure, CPLE_UserInterrupt, "User Terminated");
     443           0 :             err = CE_Failure;
     444           0 :             goto end_and_cleanup;
     445             :         }
     446             : 
     447        2124 :         err = GDALRasterIO(hRed, GF_Read, 0, iLine, nXSize, 1, pabyRedLine,
     448             :                            nXSize, 1, GDT_Byte, 0, 0);
     449        2124 :         if (err == CE_None)
     450        2124 :             err = GDALRasterIO(hGreen, GF_Read, 0, iLine, nXSize, 1,
     451             :                                pabyGreenLine, nXSize, 1, GDT_Byte, 0, 0);
     452        2124 :         if (err == CE_None)
     453        2124 :             err = GDALRasterIO(hBlue, GF_Read, 0, iLine, nXSize, 1,
     454             :                                pabyBlueLine, nXSize, 1, GDT_Byte, 0, 0);
     455        2124 :         if (err != CE_None)
     456           0 :             goto end_and_cleanup;
     457             : 
     458      479268 :         for (int iPixel = 0; iPixel < nXSize; iPixel++)
     459             :         {
     460      477144 :             const int nRed = pabyRedLine[iPixel] >> nColorShift;
     461      477144 :             const int nGreen = pabyGreenLine[iPixel] >> nColorShift;
     462      477144 :             const int nBlue = pabyBlueLine[iPixel] >> nColorShift;
     463             : 
     464      477144 :             usedboxes->rmin = std::min(usedboxes->rmin, nRed);
     465      477144 :             usedboxes->gmin = std::min(usedboxes->gmin, nGreen);
     466      477144 :             usedboxes->bmin = std::min(usedboxes->bmin, nBlue);
     467      477144 :             usedboxes->rmax = std::max(usedboxes->rmax, nRed);
     468      477144 :             usedboxes->gmax = std::max(usedboxes->gmax, nGreen);
     469      477144 :             usedboxes->bmax = std::max(usedboxes->bmax, nBlue);
     470             : 
     471             :             bool bFirstOccurrence;
     472      477144 :             if (psHashHistogram)
     473             :             {
     474      389644 :                 int *pnColor = FindAndInsertColorCount(
     475      389644 :                     psHashHistogram, MAKE_COLOR_CODE(nRed, nGreen, nBlue));
     476      389644 :                 bFirstOccurrence = (*pnColor == 0);
     477      389644 :                 (*pnColor)++;
     478             :             }
     479             :             else
     480             :             {
     481             :                 T *pnColor =
     482       87500 :                     HISTOGRAM(histogram, nCLevels, nRed, nGreen, nBlue);
     483       87500 :                 bFirstOccurrence = (*pnColor == 0);
     484       87500 :                 (*pnColor)++;
     485             :             }
     486      477144 :             if (bFirstOccurrence)
     487             :             {
     488       78587 :                 if (nColorShift == 0 && nColorCounter < nColors)
     489             :                 {
     490        2305 :                     anRed[nColorCounter] = static_cast<GByte>(nRed);
     491        2305 :                     anGreen[nColorCounter] = static_cast<GByte>(nGreen);
     492        2305 :                     anBlue[nColorCounter] = static_cast<GByte>(nBlue);
     493             :                 }
     494       78587 :                 nColorCounter++;
     495             :             }
     496             :         }
     497             :     }
     498             : 
     499          13 :     if (!pfnProgress(1.0, "Generating Histogram", pProgressArg))
     500             :     {
     501           0 :         CPLError(CE_Failure, CPLE_UserInterrupt, "User Terminated");
     502           0 :         err = CE_Failure;
     503           0 :         goto end_and_cleanup;
     504             :     }
     505             : 
     506          13 :     if (nColorShift == 0 && nColorCounter <= nColors)
     507             :     {
     508             : #if DEBUG_VERBOSE
     509             :         CPLDebug("MEDIAN_CUT", "%d colors found <= %d", nColorCounter, nColors);
     510             : #endif
     511           2 :         for (int iColor = 0; iColor < nColorCounter; iColor++)
     512             :         {
     513           1 :             const GDALColorEntry sEntry = {static_cast<GByte>(anRed[iColor]),
     514           1 :                                            static_cast<GByte>(anGreen[iColor]),
     515           1 :                                            static_cast<GByte>(anBlue[iColor]),
     516             :                                            255};
     517           1 :             GDALSetColorEntry(hColorTable, iColor, &sEntry);
     518             :         }
     519           1 :         goto end_and_cleanup;
     520             :     }
     521             : 
     522             :     /* ==================================================================== */
     523             :     /*      STEP 3: continually subdivide boxes until no more free          */
     524             :     /*      boxes remain or until all colors assigned.                      */
     525             :     /* ==================================================================== */
     526        2584 :     while (freeboxes != nullptr)
     527             :     {
     528        2572 :         auto ptr = largest_box(usedboxes);
     529        2572 :         if (ptr != nullptr)
     530        2572 :             splitbox(ptr, histogram, psHashHistogram, nCLevels, &freeboxes,
     531             :                      &usedboxes, pabyRedBand, pabyGreenBand, pabyBlueBand,
     532             :                      nPixels);
     533             :         else
     534           0 :             freeboxes = nullptr;
     535             :     }
     536             : 
     537             :     /* ==================================================================== */
     538             :     /*      STEP 4: assign colors to all boxes                              */
     539             :     /* ==================================================================== */
     540             :     {
     541          12 :         Colorbox *ptr = usedboxes;
     542        2596 :         for (int i = 0; ptr != nullptr; ++i, ptr = ptr->next)
     543             :         {
     544        2584 :             const GDALColorEntry sEntry = {
     545        2584 :                 static_cast<GByte>(((ptr->rmin + ptr->rmax) << nColorShift) /
     546             :                                    2),
     547        2584 :                 static_cast<GByte>(((ptr->gmin + ptr->gmax) << nColorShift) /
     548             :                                    2),
     549        2584 :                 static_cast<GByte>(((ptr->bmin + ptr->bmax) << nColorShift) /
     550             :                                    2),
     551             :                 255};
     552        2584 :             GDALSetColorEntry(hColorTable, i, &sEntry);
     553             :         }
     554             :     }
     555             : 
     556          12 : end_and_cleanup:
     557          13 :     CPLFree(pabyRedLine);
     558          13 :     CPLFree(pabyGreenLine);
     559          13 :     CPLFree(pabyBlueLine);
     560             : 
     561             :     // We're done with the boxes now.
     562          13 :     CPLFree(box_list);
     563          13 :     freeboxes = nullptr;
     564          13 :     usedboxes = nullptr;
     565             : 
     566          13 :     if (panHistogram == nullptr)
     567           3 :         CPLFree(histogram);
     568             : 
     569          13 :     return err;
     570             : }
     571             : 
     572             : /************************************************************************/
     573             : /*                            largest_box()                             */
     574             : /************************************************************************/
     575             : 
     576        2572 : static Colorbox *largest_box(Colorbox *usedboxes)
     577             : {
     578        2572 :     Colorbox *b = nullptr;
     579             : 
     580      329120 :     for (Colorbox *p = usedboxes; p != nullptr; p = p->next)
     581             :     {
     582      326548 :         if ((p->rmax > p->rmin || p->gmax > p->gmin || p->bmax > p->bmin) &&
     583      299299 :             (b == nullptr || p->total > b->total))
     584             :         {
     585       39803 :             b = p;
     586             :         }
     587             :     }
     588        2572 :     return b;
     589             : }
     590             : 
     591         743 : static void shrinkboxFromBand(Colorbox *ptr, const GByte *pabyRedBand,
     592             :                               const GByte *pabyGreenBand,
     593             :                               const GByte *pabyBlueBand, GUIntBig nPixels)
     594             : {
     595         743 :     int rmin_new = 255;
     596         743 :     int rmax_new = 0;
     597         743 :     int gmin_new = 255;
     598         743 :     int gmax_new = 0;
     599         743 :     int bmin_new = 255;
     600         743 :     int bmax_new = 0;
     601    35470300 :     for (GUIntBig i = 0; i < nPixels; i++)
     602             :     {
     603    35469500 :         const int iR = pabyRedBand[i];
     604    35469500 :         const int iG = pabyGreenBand[i];
     605    35469500 :         const int iB = pabyBlueBand[i];
     606    35469500 :         if (iR >= ptr->rmin && iR <= ptr->rmax && iG >= ptr->gmin &&
     607     5781920 :             iG <= ptr->gmax && iB >= ptr->bmin && iB <= ptr->bmax)
     608             :         {
     609     2154140 :             if (iR < rmin_new)
     610        3143 :                 rmin_new = iR;
     611     2154140 :             if (iR > rmax_new)
     612        4274 :                 rmax_new = iR;
     613     2154140 :             if (iG < gmin_new)
     614        3392 :                 gmin_new = iG;
     615     2154140 :             if (iG > gmax_new)
     616        4635 :                 gmax_new = iG;
     617     2154140 :             if (iB < bmin_new)
     618        4455 :                 bmin_new = iB;
     619     2154140 :             if (iB > bmax_new)
     620        3314 :                 bmax_new = iB;
     621             :         }
     622             :     }
     623             : 
     624         743 :     CPLAssert(rmin_new >= ptr->rmin && rmin_new <= rmax_new &&
     625             :               rmax_new <= ptr->rmax);
     626         743 :     CPLAssert(gmin_new >= ptr->gmin && gmin_new <= gmax_new &&
     627             :               gmax_new <= ptr->gmax);
     628         743 :     CPLAssert(bmin_new >= ptr->bmin && bmin_new <= bmax_new &&
     629             :               bmax_new <= ptr->bmax);
     630         743 :     ptr->rmin = rmin_new;
     631         743 :     ptr->rmax = rmax_new;
     632         743 :     ptr->gmin = gmin_new;
     633         743 :     ptr->gmax = gmax_new;
     634         743 :     ptr->bmin = bmin_new;
     635         743 :     ptr->bmax = bmax_new;
     636         743 : }
     637             : 
     638        3415 : static void shrinkboxFromHashHistogram(Colorbox *box,
     639             :                                        const HashHistogram *psHashHistogram)
     640             : {
     641        3415 :     if (box->rmax > box->rmin)
     642             :     {
     643        4917 :         for (int ir = box->rmin; ir <= box->rmax; ++ir)
     644             :         {
     645       54728 :             for (int ig = box->gmin; ig <= box->gmax; ++ig)
     646             :             {
     647      770529 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
     648             :                 {
     649      720718 :                     if (FindColorCount(psHashHistogram,
     650     1441440 :                                        MAKE_COLOR_CODE(ir, ig, ib)) != 0)
     651             :                     {
     652        3313 :                         box->rmin = ir;
     653        3313 :                         goto have_rmin;
     654             :                     }
     655             :                 }
     656             :             }
     657             :         }
     658             :     }
     659         102 : have_rmin:
     660        3415 :     if (box->rmax > box->rmin)
     661             :     {
     662        6050 :         for (int ir = box->rmax; ir >= box->rmin; --ir)
     663             :         {
     664       65302 :             for (int ig = box->gmin; ig <= box->gmax; ++ig)
     665             :             {
     666     1301360 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
     667             :                 {
     668     1242110 :                     if (FindColorCount(psHashHistogram,
     669     2484210 :                                        MAKE_COLOR_CODE(ir, ig, ib)) != 0)
     670             :                     {
     671        3310 :                         box->rmax = ir;
     672        3310 :                         goto have_rmax;
     673             :                     }
     674             :                 }
     675             :             }
     676             :         }
     677             :     }
     678             : 
     679         105 : have_rmax:
     680        3415 :     if (box->gmax > box->gmin)
     681             :     {
     682        7027 :         for (int ig = box->gmin; ig <= box->gmax; ++ig)
     683             :         {
     684       71453 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
     685             :             {
     686     1233930 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
     687             :                 {
     688     1169500 :                     if (FindColorCount(psHashHistogram,
     689     2339010 :                                        MAKE_COLOR_CODE(ir, ig, ib)) != 0)
     690             :                     {
     691        3252 :                         box->gmin = ig;
     692        3252 :                         goto have_gmin;
     693             :                     }
     694             :                 }
     695             :             }
     696             :         }
     697             :     }
     698             : 
     699         163 : have_gmin:
     700        3415 :     if (box->gmax > box->gmin)
     701             :     {
     702        8709 :         for (int ig = box->gmax; ig >= box->gmin; --ig)
     703             :         {
     704       99428 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
     705             :             {
     706       93961 :                 int ib = box->bmin;
     707     1785560 :                 for (; ib <= box->bmax; ++ib)
     708             :                 {
     709     1694840 :                     if (FindColorCount(psHashHistogram,
     710     3389690 :                                        MAKE_COLOR_CODE(ir, ig, ib)) != 0)
     711             :                     {
     712        3242 :                         box->gmax = ig;
     713        3242 :                         goto have_gmax;
     714             :                     }
     715             :                 }
     716             :             }
     717             :         }
     718             :     }
     719             : 
     720         173 : have_gmax:
     721        3415 :     if (box->bmax > box->bmin)
     722             :     {
     723        5004 :         for (int ib = box->bmin; ib <= box->bmax; ++ib)
     724             :         {
     725       36216 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
     726             :             {
     727      572282 :                 for (int ig = box->gmin; ig <= box->gmax; ++ig)
     728             :                 {
     729      541070 :                     if (FindColorCount(psHashHistogram,
     730     1082140 :                                        MAKE_COLOR_CODE(ir, ig, ib)) != 0)
     731             :                     {
     732        3312 :                         box->bmin = ib;
     733        3312 :                         goto have_bmin;
     734             :                     }
     735             :                 }
     736             :             }
     737             :         }
     738             :     }
     739             : 
     740         103 : have_bmin:
     741        3415 :     if (box->bmax > box->bmin)
     742             :     {
     743        6196 :         for (int ib = box->bmax; ib >= box->bmin; --ib)
     744             :         {
     745       48883 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
     746             :             {
     747      675385 :                 for (int ig = box->gmin; ig <= box->gmax; ++ig)
     748             :                 {
     749      632698 :                     if (FindColorCount(psHashHistogram,
     750     1265400 :                                        MAKE_COLOR_CODE(ir, ig, ib)) != 0)
     751             :                     {
     752        3298 :                         box->bmax = ib;
     753        3298 :                         goto have_bmax;
     754             :                     }
     755             :                 }
     756             :             }
     757             :         }
     758             :     }
     759             : 
     760         117 : have_bmax:;
     761        3415 : }
     762             : 
     763             : /************************************************************************/
     764             : /*                              splitbox()                              */
     765             : /************************************************************************/
     766             : template <class T>
     767        2572 : static void splitbox(Colorbox *ptr, const T *histogram,
     768             :                      const HashHistogram *psHashHistogram, int nCLevels,
     769             :                      Colorbox **pfreeboxes, Colorbox **pusedboxes,
     770             :                      GByte *pabyRedBand, GByte *pabyGreenBand,
     771             :                      GByte *pabyBlueBand, T nPixels)
     772             : {
     773        2572 :     T hist2[256] = {};
     774        2572 :     int first = 0;
     775        2572 :     int last = 0;
     776             : 
     777             :     enum
     778             :     {
     779             :         RED,
     780             :         GREEN,
     781             :         BLUE
     782             :     } axis;
     783             : 
     784             :     // See which axis is the largest, do a histogram along that axis.  Split at
     785             :     // median point.  Contract both new boxes to fit points and return.
     786             :     {
     787        2572 :         int i = ptr->rmax - ptr->rmin;
     788        2572 :         if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
     789         858 :             axis = RED;
     790        1714 :         else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
     791         785 :             axis = GREEN;
     792             :         else
     793         929 :             axis = BLUE;
     794             :     }
     795             :     // Get histogram along longest axis.
     796        2572 :     const GUInt32 nIters = (ptr->rmax - ptr->rmin + 1) *
     797        2572 :                            (ptr->gmax - ptr->gmin + 1) *
     798        2572 :                            (ptr->bmax - ptr->bmin + 1);
     799             : 
     800        2572 :     switch (axis)
     801             :     {
     802         858 :         case RED:
     803             :         {
     804         858 :             if (nPixels != 0 && nIters > nPixels)
     805             :             {
     806         229 :                 const int rmin = ptr->rmin;
     807         229 :                 const int rmax = ptr->rmax;
     808         229 :                 const int gmin = ptr->gmin;
     809         229 :                 const int gmax = ptr->gmax;
     810         229 :                 const int bmin = ptr->bmin;
     811         229 :                 const int bmax = ptr->bmax;
     812    11014200 :                 for (T iPixel = 0; iPixel < nPixels; iPixel++)
     813             :                 {
     814    11014000 :                     int iR = pabyRedBand[iPixel];
     815    11014000 :                     int iG = pabyGreenBand[iPixel];
     816    11014000 :                     int iB = pabyBlueBand[iPixel];
     817    11014000 :                     if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
     818     1379080 :                         iB >= bmin && iB <= bmax)
     819             :                     {
     820      893690 :                         hist2[iR]++;
     821             :                     }
     822         229 :                 }
     823             :             }
     824         629 :             else if (psHashHistogram)
     825             :             {
     826         472 :                 T *histp = &hist2[ptr->rmin];
     827       11243 :                 for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
     828             :                 {
     829       10771 :                     *histp = 0;
     830      226089 :                     for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
     831             :                     {
     832     3704290 :                         for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
     833             :                         {
     834     3488980 :                             *histp += FindColorCount(
     835     3488980 :                                 psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
     836             :                         }
     837             :                     }
     838       10771 :                     histp++;
     839             :                 }
     840             :             }
     841             :             else
     842             :             {
     843         157 :                 T *histp = &hist2[ptr->rmin];
     844        2449 :                 for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
     845             :                 {
     846        2292 :                     *histp = 0;
     847       51509 :                     for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
     848             :                     {
     849             :                         const T *iptr =
     850       49217 :                             HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
     851     1135910 :                         for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
     852     1086700 :                             *histp += *iptr++;
     853             :                     }
     854        2292 :                     histp++;
     855             :                 }
     856             :             }
     857         858 :             first = ptr->rmin;
     858         858 :             last = ptr->rmax;
     859         858 :             break;
     860             :         }
     861         785 :         case GREEN:
     862             :         {
     863         785 :             if (nPixels != 0 && nIters > nPixels)
     864             :             {
     865         153 :                 const int rmin = ptr->rmin;
     866         153 :                 const int rmax = ptr->rmax;
     867         153 :                 const int gmin = ptr->gmin;
     868         153 :                 const int gmax = ptr->gmax;
     869         153 :                 const int bmin = ptr->bmin;
     870         153 :                 const int bmax = ptr->bmax;
     871     7442490 :                 for (T iPixel = 0; iPixel < nPixels; iPixel++)
     872             :                 {
     873     7442340 :                     const int iR = pabyRedBand[iPixel];
     874     7442340 :                     const int iG = pabyGreenBand[iPixel];
     875     7442340 :                     const int iB = pabyBlueBand[iPixel];
     876     7442340 :                     if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
     877     1575360 :                         iB >= bmin && iB <= bmax)
     878             :                     {
     879      948439 :                         hist2[iG]++;
     880             :                     }
     881         153 :                 }
     882             :             }
     883         632 :             else if (psHashHistogram)
     884             :             {
     885         476 :                 T *histp = &hist2[ptr->gmin];
     886       12600 :                 for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
     887             :                 {
     888       12124 :                     *histp = 0;
     889      183557 :                     for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
     890             :                     {
     891     3551710 :                         for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
     892             :                         {
     893     3380280 :                             *histp += FindColorCount(
     894     3380280 :                                 psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
     895             :                         }
     896             :                     }
     897       12124 :                     histp++;
     898             :                 }
     899             :             }
     900             :             else
     901             :             {
     902         156 :                 T *histp = &hist2[ptr->gmin];
     903        2610 :                 for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
     904             :                 {
     905        2454 :                     *histp = 0;
     906       40620 :                     for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
     907             :                     {
     908             :                         const T *iptr =
     909       38166 :                             HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
     910      851066 :                         for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
     911      812900 :                             *histp += *iptr++;
     912             :                     }
     913        2454 :                     histp++;
     914             :                 }
     915             :             }
     916         785 :             first = ptr->gmin;
     917         785 :             last = ptr->gmax;
     918         785 :             break;
     919             :         }
     920         929 :         case BLUE:
     921             :         {
     922         929 :             if (nPixels != 0 && nIters > nPixels)
     923             :             {
     924         199 :                 const int rmin = ptr->rmin;
     925         199 :                 const int rmax = ptr->rmax;
     926         199 :                 const int gmin = ptr->gmin;
     927         199 :                 const int gmax = ptr->gmax;
     928         199 :                 const int bmin = ptr->bmin;
     929         199 :                 const int bmax = ptr->bmax;
     930     9614150 :                 for (T iPixel = 0; iPixel < nPixels; iPixel++)
     931             :                 {
     932     9613950 :                     const int iR = pabyRedBand[iPixel];
     933     9613950 :                     const int iG = pabyGreenBand[iPixel];
     934     9613950 :                     const int iB = pabyBlueBand[iPixel];
     935     9613950 :                     if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
     936      832271 :                         iB >= bmin && iB <= bmax)
     937             :                     {
     938      649243 :                         hist2[iB]++;
     939             :                     }
     940         199 :                 }
     941             :             }
     942         730 :             else if (psHashHistogram)
     943             :             {
     944         573 :                 T *histp = &hist2[ptr->bmin];
     945       16540 :                 for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
     946             :                 {
     947       15967 :                     *histp = 0;
     948      252363 :                     for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
     949             :                     {
     950     4915820 :                         for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
     951             :                         {
     952     4679430 :                             *histp += FindColorCount(
     953     4679430 :                                 psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
     954             :                         }
     955             :                     }
     956       15967 :                     histp++;
     957             :                 }
     958             :             }
     959             :             else
     960             :             {
     961         157 :                 T *histp = &hist2[ptr->bmin];
     962        3487 :                 for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
     963             :                 {
     964        3330 :                     *histp = 0;
     965       57185 :                     for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
     966             :                     {
     967             :                         const T *iptr =
     968       53855 :                             HISTOGRAM(histogram, nCLevels, ir, ptr->gmin, ib);
     969     1396220 :                         for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
     970             :                         {
     971     1342360 :                             *histp += *iptr;
     972     1342360 :                             iptr += nCLevels;
     973             :                         }
     974             :                     }
     975        3330 :                     histp++;
     976             :                 }
     977             :             }
     978         929 :             first = ptr->bmin;
     979         929 :             last = ptr->bmax;
     980         929 :             break;
     981             :         }
     982             :     }
     983             :     // Find median point.
     984        2572 :     T *histp = &hist2[first];
     985        2572 :     int i = first;  // TODO(schwehr): Rename i.
     986             :     {
     987        2572 :         T sum = 0;
     988        2572 :         T sum2 = static_cast<T>(ptr->total / 2);
     989       42342 :         for (; i <= last && (sum += *histp++) < sum2; ++i)
     990             :         {
     991             :         }
     992             :     }
     993        2572 :     if (i == first)
     994         191 :         i++;
     995             : 
     996             :     // Create new box, re-allocate points.
     997        2572 :     Colorbox *new_cb = *pfreeboxes;
     998        2572 :     *pfreeboxes = new_cb->next;
     999        2572 :     if (*pfreeboxes)
    1000        2560 :         (*pfreeboxes)->prev = nullptr;
    1001        2572 :     if (*pusedboxes)
    1002        2572 :         (*pusedboxes)->prev = new_cb;
    1003        2572 :     new_cb->next = *pusedboxes;
    1004        2572 :     *pusedboxes = new_cb;
    1005             : 
    1006        2572 :     histp = &hist2[first];
    1007             :     {
    1008        2572 :         T sum1 = 0;
    1009       42533 :         for (int j = first; j < i; j++)
    1010       39961 :             sum1 += *histp++;
    1011        2572 :         T sum2 = 0;
    1012       65856 :         for (int j = i; j <= last; j++)
    1013       63284 :             sum2 += *histp++;
    1014        2572 :         new_cb->total = sum1;
    1015        2572 :         ptr->total = sum2;
    1016             :     }
    1017        2572 :     new_cb->rmin = ptr->rmin;
    1018        2572 :     new_cb->rmax = ptr->rmax;
    1019        2572 :     new_cb->gmin = ptr->gmin;
    1020        2572 :     new_cb->gmax = ptr->gmax;
    1021        2572 :     new_cb->bmin = ptr->bmin;
    1022        2572 :     new_cb->bmax = ptr->bmax;
    1023        2572 :     switch (axis)
    1024             :     {
    1025         858 :         case RED:
    1026         858 :             new_cb->rmax = i - 1;
    1027         858 :             ptr->rmin = i;
    1028         858 :             break;
    1029         785 :         case GREEN:
    1030         785 :             new_cb->gmax = i - 1;
    1031         785 :             ptr->gmin = i;
    1032         785 :             break;
    1033         929 :         case BLUE:
    1034         929 :             new_cb->bmax = i - 1;
    1035         929 :             ptr->bmin = i;
    1036         929 :             break;
    1037             :     }
    1038             : 
    1039        2572 :     if (nPixels != 0 &&
    1040        2295 :         static_cast<T>(new_cb->rmax - new_cb->rmin + 1) *
    1041        2295 :                 static_cast<T>(new_cb->gmax - new_cb->gmin + 1) *
    1042        2295 :                 static_cast<T>(new_cb->bmax - new_cb->bmin + 1) >
    1043             :             nPixels)
    1044             :     {
    1045         325 :         shrinkboxFromBand(new_cb, pabyRedBand, pabyGreenBand, pabyBlueBand,
    1046             :                           nPixels);
    1047             :     }
    1048        2247 :     else if (psHashHistogram != nullptr)
    1049             :     {
    1050        1748 :         shrinkboxFromHashHistogram(new_cb, psHashHistogram);
    1051             :     }
    1052             :     else
    1053             :     {
    1054         499 :         shrinkbox(new_cb, histogram, nCLevels);
    1055             :     }
    1056             : 
    1057        2572 :     if (nPixels != 0 && static_cast<T>(ptr->rmax - ptr->rmin + 1) *
    1058        2295 :                                 static_cast<T>(ptr->gmax - ptr->gmin + 1) *
    1059        2295 :                                 static_cast<T>(ptr->bmax - ptr->bmin + 1) >
    1060             :                             nPixels)
    1061             :     {
    1062         418 :         shrinkboxFromBand(ptr, pabyRedBand, pabyGreenBand, pabyBlueBand,
    1063             :                           nPixels);
    1064             :     }
    1065        2154 :     else if (psHashHistogram != nullptr)
    1066             :     {
    1067        1667 :         shrinkboxFromHashHistogram(ptr, psHashHistogram);
    1068             :     }
    1069             :     else
    1070             :     {
    1071         487 :         shrinkbox(ptr, histogram, nCLevels);
    1072             :     }
    1073        2572 : }
    1074             : 
    1075             : /************************************************************************/
    1076             : /*                             shrinkbox()                              */
    1077             : /************************************************************************/
    1078             : template <class T>
    1079         986 : static void shrinkbox(Colorbox *box, const T *histogram, int nCLevels)
    1080             : {
    1081         986 :     if (box->rmax > box->rmin)
    1082             :     {
    1083         972 :         for (int ir = box->rmin; ir <= box->rmax; ++ir)
    1084             :         {
    1085        7471 :             for (int ig = box->gmin; ig <= box->gmax; ++ig)
    1086             :             {
    1087             :                 const T *histp =
    1088        7202 :                     HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
    1089      162688 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
    1090             :                 {
    1091      156189 :                     if (*histp++ != 0)
    1092             :                     {
    1093         703 :                         box->rmin = ir;
    1094         703 :                         goto have_rmin;
    1095             :                     }
    1096             :                 }
    1097             :             }
    1098             :         }
    1099             :     }
    1100         283 : have_rmin:
    1101         986 :     if (box->rmax > box->rmin)
    1102             :     {
    1103        1384 :         for (int ir = box->rmax; ir >= box->rmin; --ir)
    1104             :         {
    1105       10381 :             for (int ig = box->gmin; ig <= box->gmax; ++ig)
    1106             :             {
    1107             :                 const T *histp =
    1108        9691 :                     HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
    1109      261708 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
    1110             :                 {
    1111      252711 :                     if (*histp++ != 0)
    1112             :                     {
    1113         694 :                         box->rmax = ir;
    1114         694 :                         goto have_rmax;
    1115             :                     }
    1116             :                 }
    1117             :             }
    1118             :         }
    1119             :     }
    1120             : 
    1121         292 : have_rmax:
    1122         986 :     if (box->gmax > box->gmin)
    1123             :     {
    1124        1173 :         for (int ig = box->gmin; ig <= box->gmax; ++ig)
    1125             :         {
    1126       10451 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
    1127             :             {
    1128             :                 const T *histp =
    1129        9963 :                     HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
    1130      251615 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
    1131             :                 {
    1132      242337 :                     if (*histp++ != 0)
    1133             :                     {
    1134         685 :                         box->gmin = ig;
    1135         685 :                         goto have_gmin;
    1136             :                     }
    1137             :                 }
    1138             :             }
    1139             :         }
    1140             :     }
    1141             : 
    1142         301 : have_gmin:
    1143         986 :     if (box->gmax > box->gmin)
    1144             :     {
    1145        1598 :         for (int ig = box->gmax; ig >= box->gmin; --ig)
    1146             :         {
    1147       17164 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
    1148             :             {
    1149             :                 const T *histp =
    1150       16244 :                     HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
    1151      353961 :                 for (int ib = box->bmin; ib <= box->bmax; ++ib)
    1152             :                 {
    1153      338395 :                     if (*histp++ != 0)
    1154             :                     {
    1155         678 :                         box->gmax = ig;
    1156         678 :                         goto have_gmax;
    1157             :                     }
    1158             :                 }
    1159             :             }
    1160             :         }
    1161             :     }
    1162             : 
    1163         308 : have_gmax:
    1164         986 :     if (box->bmax > box->bmin)
    1165             :     {
    1166         939 :         for (int ib = box->bmin; ib <= box->bmax; ++ib)
    1167             :         {
    1168        5880 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
    1169             :             {
    1170             :                 const T *histp =
    1171        5642 :                     HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
    1172      102928 :                 for (int ig = box->gmin; ig <= box->gmax; ++ig)
    1173             :                 {
    1174       97987 :                     if (*histp != 0)
    1175             :                     {
    1176         701 :                         box->bmin = ib;
    1177         701 :                         goto have_bmin;
    1178             :                     }
    1179       97286 :                     histp += nCLevels;
    1180             :                 }
    1181             :             }
    1182             :         }
    1183             :     }
    1184             : 
    1185         285 : have_bmin:
    1186         986 :     if (box->bmax > box->bmin)
    1187             :     {
    1188        1341 :         for (int ib = box->bmax; ib >= box->bmin; --ib)
    1189             :         {
    1190       10742 :             for (int ir = box->rmin; ir <= box->rmax; ++ir)
    1191             :             {
    1192             :                 const T *histp =
    1193       10083 :                     HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
    1194      112931 :                 for (int ig = box->gmin; ig <= box->gmax; ++ig)
    1195             :                 {
    1196      103530 :                     if (*histp != 0)
    1197             :                     {
    1198         682 :                         box->bmax = ib;
    1199         682 :                         goto have_bmax;
    1200             :                     }
    1201      102848 :                     histp += nCLevels;
    1202             :                 }
    1203             :             }
    1204             :         }
    1205             :     }
    1206             : 
    1207         304 : have_bmax:;
    1208         986 : }
    1209             : 
    1210             : // Explicitly instantiate template functions.
    1211             : template int GDALComputeMedianCutPCTInternal<GUInt32>(
    1212             :     GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
    1213             :     GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
    1214             :     int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
    1215             :     GUInt32 *panHistogram, GDALColorTableH hColorTable,
    1216             :     GDALProgressFunc pfnProgress, void *pProgressArg);
    1217             : 
    1218             : template int GDALComputeMedianCutPCTInternal<GUIntBig>(
    1219             :     GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
    1220             :     GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
    1221             :     int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
    1222             :     GUIntBig *panHistogram, GDALColorTableH hColorTable,
    1223             :     GDALProgressFunc pfnProgress, void *pProgressArg);

Generated by: LCOV version 1.14