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 4128224 : template <typename T> static T *HISTOGRAM(T *h, int n, int r, int g, int b)
42 : {
43 4128224 : const int index = (r * n + g) * n + b;
44 4128224 : 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 : * Since GDAL 3.13, source nodata values or mask band will be taken into account
116 : * to determine which pixels are valid.
117 : *
118 : * @param hRed Red input band.
119 : * @param hGreen Green input band.
120 : * @param hBlue Blue input band.
121 : * @param pfnIncludePixel function used to test which pixels should be included
122 : * in the analysis. At this time this argument is ignored.
123 : * This should normally be NULL.
124 : * @param nColors the desired number of colors to be returned (2-256).
125 : * @param hColorTable the colors will be returned in this color table object.
126 : * @param pfnProgress callback for reporting algorithm progress matching the
127 : * GDALProgressFunc() semantics. May be NULL.
128 : * @param pProgressArg callback argument passed to pfnProgress.
129 : *
130 : * @return returns CE_None on success or CE_Failure if an error occurs.
131 : */
132 :
133 4 : extern "C" int CPL_STDCALL GDALComputeMedianCutPCT(
134 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
135 : int (*pfnIncludePixel)(int, int, void *), int nColors,
136 : GDALColorTableH hColorTable, GDALProgressFunc pfnProgress,
137 : void *pProgressArg)
138 :
139 : {
140 4 : VALIDATE_POINTER1(hRed, "GDALComputeMedianCutPCT", CE_Failure);
141 4 : const int nXSize = GDALGetRasterBandXSize(hRed);
142 4 : const int nYSize = GDALGetRasterBandYSize(hRed);
143 4 : if (nYSize == 0)
144 0 : return CE_Failure;
145 4 : if (static_cast<GUInt32>(nXSize) <
146 4 : std::numeric_limits<GUInt32>::max() / static_cast<GUInt32>(nYSize))
147 : {
148 4 : return GDALComputeMedianCutPCTInternal(
149 : hRed, hGreen, hBlue, nullptr, nullptr, nullptr, pfnIncludePixel,
150 : nColors, 5, static_cast<GUInt32 *>(nullptr), hColorTable,
151 4 : pfnProgress, pProgressArg);
152 : }
153 : else
154 : {
155 0 : return GDALComputeMedianCutPCTInternal(
156 : hRed, hGreen, hBlue, nullptr, nullptr, nullptr, pfnIncludePixel,
157 : nColors, 5, static_cast<GUIntBig *>(nullptr), hColorTable,
158 0 : pfnProgress, pProgressArg);
159 : }
160 : }
161 :
162 : #ifndef IsColorCodeSet_defined
163 : #define IsColorCodeSet_defined
164 :
165 18835800 : static inline bool IsColorCodeSet(GUInt32 nColorCode)
166 : {
167 18835800 : return (nColorCode >> 31) == 0;
168 : }
169 : #endif
170 :
171 17549600 : static inline int FindColorCount(const HashHistogram *psHashHistogram,
172 : GUInt32 nColorCode)
173 : {
174 :
175 17549600 : GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
176 : while (true)
177 : {
178 17550200 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode))
179 : {
180 16219100 : return 0;
181 : }
182 1331150 : if (psHashHistogram[nIdx].nColorCode == nColorCode)
183 : {
184 141464 : return psHashHistogram[nIdx].nCount;
185 : }
186 1189680 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode2))
187 : {
188 1151920 : return 0;
189 : }
190 37761 : if (psHashHistogram[nIdx].nColorCode2 == nColorCode)
191 : {
192 4030 : return psHashHistogram[nIdx].nCount2;
193 : }
194 33731 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode3))
195 : {
196 33078 : return 0;
197 : }
198 653 : if (psHashHistogram[nIdx].nColorCode3 == nColorCode)
199 : {
200 59 : return psHashHistogram[nIdx].nCount3;
201 : }
202 :
203 0 : do
204 : {
205 594 : nIdx += 257;
206 594 : if (nIdx >= PRIME_FOR_65536)
207 0 : nIdx -= PRIME_FOR_65536;
208 594 : } while (IsColorCodeSet(psHashHistogram[nIdx].nColorCode) &&
209 260 : psHashHistogram[nIdx].nColorCode != nColorCode &&
210 130 : IsColorCodeSet(psHashHistogram[nIdx].nColorCode2) &&
211 0 : psHashHistogram[nIdx].nColorCode2 != nColorCode &&
212 724 : IsColorCodeSet(psHashHistogram[nIdx].nColorCode3) &&
213 0 : psHashHistogram[nIdx].nColorCode3 != nColorCode);
214 : }
215 : }
216 :
217 389644 : static inline int *FindAndInsertColorCount(HashHistogram *psHashHistogram,
218 : GUInt32 nColorCode)
219 : {
220 389644 : GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
221 : while (true)
222 : {
223 389644 : if (psHashHistogram[nIdx].nColorCode == nColorCode)
224 : {
225 329957 : return &(psHashHistogram[nIdx].nCount);
226 : }
227 59687 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode))
228 : {
229 56306 : psHashHistogram[nIdx].nColorCode = nColorCode;
230 56306 : psHashHistogram[nIdx].nCount = 0;
231 56306 : return &(psHashHistogram[nIdx].nCount);
232 : }
233 3381 : if (psHashHistogram[nIdx].nColorCode2 == nColorCode)
234 : {
235 1661 : return &(psHashHistogram[nIdx].nCount2);
236 : }
237 1720 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode2))
238 : {
239 1687 : psHashHistogram[nIdx].nColorCode2 = nColorCode;
240 1687 : psHashHistogram[nIdx].nCount2 = 0;
241 1687 : return &(psHashHistogram[nIdx].nCount2);
242 : }
243 33 : if (psHashHistogram[nIdx].nColorCode3 == nColorCode)
244 : {
245 5 : return &(psHashHistogram[nIdx].nCount3);
246 : }
247 28 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode3))
248 : {
249 28 : psHashHistogram[nIdx].nColorCode3 = nColorCode;
250 28 : psHashHistogram[nIdx].nCount3 = 0;
251 28 : return &(psHashHistogram[nIdx].nCount3);
252 : }
253 :
254 0 : do
255 : {
256 0 : nIdx += 257;
257 0 : if (nIdx >= PRIME_FOR_65536)
258 0 : nIdx -= PRIME_FOR_65536;
259 0 : } while (IsColorCodeSet(psHashHistogram[nIdx].nColorCode) &&
260 0 : psHashHistogram[nIdx].nColorCode != nColorCode &&
261 0 : IsColorCodeSet(psHashHistogram[nIdx].nColorCode2) &&
262 0 : psHashHistogram[nIdx].nColorCode2 != nColorCode &&
263 0 : IsColorCodeSet(psHashHistogram[nIdx].nColorCode3) &&
264 0 : psHashHistogram[nIdx].nColorCode3 != nColorCode);
265 : }
266 : }
267 :
268 : template <class T>
269 50 : int GDALComputeMedianCutPCTInternal(
270 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
271 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
272 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
273 : T *panHistogram, // NULL, or >= size (1<<nBits)^3 * sizeof(T) bytes.
274 : GDALColorTableH hColorTable, GDALProgressFunc pfnProgress,
275 : void *pProgressArg, std::vector<T> *panPixelCountPerColorTableEntry)
276 :
277 : {
278 50 : VALIDATE_POINTER1(hRed, "GDALComputeMedianCutPCT", CE_Failure);
279 50 : VALIDATE_POINTER1(hGreen, "GDALComputeMedianCutPCT", CE_Failure);
280 50 : VALIDATE_POINTER1(hBlue, "GDALComputeMedianCutPCT", CE_Failure);
281 :
282 50 : CPLErr err = CE_None;
283 :
284 : /* -------------------------------------------------------------------- */
285 : /* Validate parameters. */
286 : /* -------------------------------------------------------------------- */
287 50 : const int nXSize = GDALGetRasterBandXSize(hRed);
288 50 : const int nYSize = GDALGetRasterBandYSize(hRed);
289 :
290 50 : if (GDALGetRasterBandXSize(hGreen) != nXSize ||
291 50 : GDALGetRasterBandYSize(hGreen) != nYSize ||
292 150 : GDALGetRasterBandXSize(hBlue) != nXSize ||
293 50 : GDALGetRasterBandYSize(hBlue) != nYSize)
294 : {
295 0 : CPLError(CE_Failure, CPLE_IllegalArg,
296 : "Green or blue band doesn't match size of red band.");
297 :
298 0 : return CE_Failure;
299 : }
300 :
301 50 : if (pfnIncludePixel != nullptr)
302 : {
303 0 : CPLError(CE_Failure, CPLE_IllegalArg,
304 : "GDALComputeMedianCutPCT() doesn't currently support "
305 : "pfnIncludePixel function.");
306 :
307 0 : return CE_Failure;
308 : }
309 :
310 50 : if (nColors <= 0)
311 : {
312 0 : CPLError(CE_Failure, CPLE_IllegalArg,
313 : "GDALComputeMedianCutPCT(): "
314 : "nColors must be strictly greater than 1.");
315 :
316 0 : return CE_Failure;
317 : }
318 :
319 50 : if (nColors > 256)
320 : {
321 0 : CPLError(CE_Failure, CPLE_IllegalArg,
322 : "GDALComputeMedianCutPCT(): "
323 : "nColors must be lesser than or equal to 256.");
324 :
325 0 : return CE_Failure;
326 : }
327 :
328 50 : if (pfnProgress == nullptr)
329 23 : pfnProgress = GDALDummyProgress;
330 :
331 50 : int nSrcNoData = -1;
332 50 : GDALRasterBandH hMaskBand = nullptr;
333 50 : int bSrcHasNoDataR = FALSE;
334 50 : const double dfSrcNoDataR = GDALGetRasterNoDataValue(hRed, &bSrcHasNoDataR);
335 50 : int bSrcHasNoDataG = FALSE;
336 : const double dfSrcNoDataG =
337 50 : GDALGetRasterNoDataValue(hGreen, &bSrcHasNoDataG);
338 50 : int bSrcHasNoDataB = FALSE;
339 : const double dfSrcNoDataB =
340 50 : GDALGetRasterNoDataValue(hBlue, &bSrcHasNoDataB);
341 50 : if (bSrcHasNoDataR && bSrcHasNoDataG && bSrcHasNoDataB &&
342 2 : dfSrcNoDataR == dfSrcNoDataG && dfSrcNoDataR == dfSrcNoDataB &&
343 2 : dfSrcNoDataR >= 0 && dfSrcNoDataR <= 255 &&
344 2 : std::round(dfSrcNoDataR) == dfSrcNoDataR)
345 : {
346 2 : nSrcNoData = static_cast<int>(dfSrcNoDataR);
347 : }
348 : else
349 : {
350 48 : const int nMaskFlags = GDALGetMaskFlags(hRed);
351 48 : if ((nMaskFlags & GMF_PER_DATASET))
352 3 : hMaskBand = GDALGetMaskBand(hRed);
353 : }
354 :
355 : /* ==================================================================== */
356 : /* STEP 1: create empty boxes. */
357 : /* ==================================================================== */
358 71 : if (static_cast<GUInt32>(nXSize) >
359 50 : cpl::NumericLimits<T>::max() / static_cast<GUInt32>(nYSize))
360 : {
361 0 : CPLError(CE_Warning, CPLE_AppDefined,
362 : "GDALComputeMedianCutPCTInternal() not called "
363 : "with large enough type");
364 : }
365 :
366 50 : T nPixels = 0;
367 11 : if (nBits == 8 && pabyRedBand != nullptr && pabyGreenBand != nullptr &&
368 61 : pabyBlueBand != nullptr &&
369 10 : static_cast<GUInt32>(nXSize) <=
370 10 : cpl::NumericLimits<T>::max() / static_cast<GUInt32>(nYSize))
371 : {
372 10 : nPixels = static_cast<T>(nXSize) * static_cast<T>(nYSize);
373 : }
374 :
375 50 : const int nCLevels = 1 << nBits;
376 50 : const int nCLevelsCube = nCLevels * nCLevels * nCLevels;
377 50 : T *histogram = nullptr;
378 50 : HashHistogram *psHashHistogram = nullptr;
379 50 : if (panHistogram)
380 : {
381 10 : if (nBits == 8 && static_cast<GUIntBig>(nXSize) * nYSize <= 65536)
382 : {
383 : // If the image is small enough, then the number of colors
384 : // will be limited and using a hashmap, rather than a full table
385 : // will be more efficient.
386 9 : histogram = nullptr;
387 9 : psHashHistogram = reinterpret_cast<HashHistogram *>(panHistogram);
388 9 : memset(psHashHistogram, 0xFF,
389 : sizeof(HashHistogram) * PRIME_FOR_65536);
390 : }
391 : else
392 : {
393 1 : histogram = panHistogram;
394 1 : memset(histogram, 0, nCLevelsCube * sizeof(T));
395 : }
396 : }
397 : else
398 : {
399 : histogram =
400 40 : static_cast<T *>(VSI_CALLOC_VERBOSE(nCLevelsCube, sizeof(T)));
401 40 : if (histogram == nullptr)
402 : {
403 0 : return CE_Failure;
404 : }
405 : }
406 : Colorbox *box_list =
407 50 : static_cast<Colorbox *>(CPLMalloc(nColors * sizeof(Colorbox)));
408 50 : Colorbox *freeboxes = box_list;
409 50 : freeboxes[0].next = &freeboxes[1];
410 50 : freeboxes[0].prev = nullptr;
411 11339 : for (int i = 1; i < nColors - 1; ++i)
412 : {
413 11289 : freeboxes[i].next = &freeboxes[i + 1];
414 11289 : freeboxes[i].prev = &freeboxes[i - 1];
415 : }
416 50 : freeboxes[nColors - 1].next = nullptr;
417 50 : freeboxes[nColors - 1].prev = &freeboxes[nColors - 2];
418 :
419 : /* ==================================================================== */
420 : /* Build histogram. */
421 : /* ==================================================================== */
422 :
423 : /* -------------------------------------------------------------------- */
424 : /* Initialize the box datastructures. */
425 : /* -------------------------------------------------------------------- */
426 50 : Colorbox *freeboxes_before = freeboxes;
427 50 : freeboxes = freeboxes_before->next;
428 50 : if (freeboxes)
429 50 : freeboxes->prev = nullptr;
430 :
431 50 : Colorbox *usedboxes = freeboxes_before;
432 50 : usedboxes->next = nullptr;
433 50 : usedboxes->rmin = 999;
434 50 : usedboxes->gmin = 999;
435 50 : usedboxes->bmin = 999;
436 50 : usedboxes->rmax = -1;
437 50 : usedboxes->gmax = -1;
438 50 : usedboxes->bmax = -1;
439 50 : usedboxes->total =
440 50 : static_cast<GUIntBig>(nXSize) * static_cast<GUIntBig>(nYSize);
441 :
442 : /* -------------------------------------------------------------------- */
443 : /* Collect histogram. */
444 : /* -------------------------------------------------------------------- */
445 :
446 : // TODO(schwehr): Move these closer to usage after removing gotos.
447 50 : const int nColorShift = 8 - nBits;
448 50 : int nColorCounter = 0;
449 50 : GByte anRed[256] = {};
450 50 : GByte anGreen[256] = {};
451 50 : GByte anBlue[256] = {};
452 :
453 50 : GByte *pabyRedLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
454 50 : GByte *pabyGreenLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
455 50 : GByte *pabyBlueLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
456 0 : std::unique_ptr<GByte, VSIFreeReleaser> pabyMask;
457 50 : if (hMaskBand)
458 3 : pabyMask.reset(static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)));
459 :
460 50 : if (pabyRedLine == nullptr || pabyGreenLine == nullptr ||
461 100 : pabyBlueLine == nullptr || (hMaskBand && !pabyMask))
462 : {
463 0 : err = CE_Failure;
464 0 : goto end_and_cleanup;
465 : }
466 :
467 8549 : for (int iLine = 0; iLine < nYSize; iLine++)
468 : {
469 8499 : if (!pfnProgress(iLine / static_cast<double>(nYSize),
470 : "Generating Histogram", pProgressArg))
471 : {
472 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User Terminated");
473 0 : err = CE_Failure;
474 0 : goto end_and_cleanup;
475 : }
476 :
477 8499 : err = GDALRasterIO(hRed, GF_Read, 0, iLine, nXSize, 1, pabyRedLine,
478 : nXSize, 1, GDT_UInt8, 0, 0);
479 8499 : if (err == CE_None)
480 8499 : err = GDALRasterIO(hGreen, GF_Read, 0, iLine, nXSize, 1,
481 : pabyGreenLine, nXSize, 1, GDT_UInt8, 0, 0);
482 8499 : if (err == CE_None)
483 8499 : err = GDALRasterIO(hBlue, GF_Read, 0, iLine, nXSize, 1,
484 : pabyBlueLine, nXSize, 1, GDT_UInt8, 0, 0);
485 8499 : if (err == CE_None && hMaskBand)
486 1687 : err = GDALRasterIO(hMaskBand, GF_Read, 0, iLine, nXSize, 1,
487 1687 : pabyMask.get(), nXSize, 1, GDT_UInt8, 0, 0);
488 8499 : if (err != CE_None)
489 0 : goto end_and_cleanup;
490 :
491 6324100 : for (int iPixel = 0; iPixel < nXSize; iPixel++)
492 : {
493 15038530 : if ((pabyRedLine[iPixel] == nSrcNoData &&
494 33554 : pabyGreenLine[iPixel] == nSrcNoData &&
495 15014790 : pabyBlueLine[iPixel] == nSrcNoData) ||
496 8665890 : (pabyMask && (pabyMask.get())[iPixel] == 0))
497 : {
498 2373738 : continue;
499 : }
500 :
501 3941870 : const int nRed = pabyRedLine[iPixel] >> nColorShift;
502 3941870 : const int nGreen = pabyGreenLine[iPixel] >> nColorShift;
503 3941870 : const int nBlue = pabyBlueLine[iPixel] >> nColorShift;
504 :
505 3941870 : usedboxes->rmin = std::min(usedboxes->rmin, nRed);
506 3941870 : usedboxes->gmin = std::min(usedboxes->gmin, nGreen);
507 3941870 : usedboxes->bmin = std::min(usedboxes->bmin, nBlue);
508 3941870 : usedboxes->rmax = std::max(usedboxes->rmax, nRed);
509 3941870 : usedboxes->gmax = std::max(usedboxes->gmax, nGreen);
510 3941870 : usedboxes->bmax = std::max(usedboxes->bmax, nBlue);
511 :
512 : bool bFirstOccurrence;
513 3941870 : if (psHashHistogram)
514 : {
515 389644 : int *pnColor = FindAndInsertColorCount(
516 389644 : psHashHistogram, MAKE_COLOR_CODE(nRed, nGreen, nBlue));
517 389644 : bFirstOccurrence = (*pnColor == 0);
518 389644 : (*pnColor)++;
519 : }
520 : else
521 : {
522 : T *pnColor =
523 3552220 : HISTOGRAM(histogram, nCLevels, nRed, nGreen, nBlue);
524 3552220 : bFirstOccurrence = (*pnColor == 0);
525 3552220 : (*pnColor)++;
526 : }
527 3941870 : if (bFirstOccurrence)
528 : {
529 111366 : if (nColorShift == 0 && nColorCounter < nColors)
530 : {
531 2314 : anRed[nColorCounter] = static_cast<GByte>(nRed);
532 2314 : anGreen[nColorCounter] = static_cast<GByte>(nGreen);
533 2314 : anBlue[nColorCounter] = static_cast<GByte>(nBlue);
534 : }
535 111366 : nColorCounter++;
536 : }
537 : }
538 : }
539 :
540 50 : if (!pfnProgress(1.0, "Generating Histogram", pProgressArg))
541 : {
542 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User Terminated");
543 0 : err = CE_Failure;
544 0 : goto end_and_cleanup;
545 : }
546 :
547 50 : if (nColorShift == 0 && nColorCounter <= nColors)
548 : {
549 : #if DEBUG_VERBOSE
550 : CPLDebug("MEDIAN_CUT", "%d colors found <= %d", nColorCounter, nColors);
551 : #endif
552 2 : if (panPixelCountPerColorTableEntry)
553 : {
554 1 : panPixelCountPerColorTableEntry->clear();
555 1 : panPixelCountPerColorTableEntry->reserve(nColors);
556 : }
557 12 : for (int iColor = 0; iColor < nColorCounter; iColor++)
558 : {
559 10 : const GDALColorEntry sEntry = {static_cast<GByte>(anRed[iColor]),
560 10 : static_cast<GByte>(anGreen[iColor]),
561 10 : static_cast<GByte>(anBlue[iColor]),
562 : 255};
563 10 : GDALSetColorEntry(hColorTable, iColor, &sEntry);
564 10 : if (panPixelCountPerColorTableEntry)
565 : {
566 9 : if (psHashHistogram)
567 : {
568 0 : int *pnColor = FindAndInsertColorCount(
569 : psHashHistogram,
570 0 : MAKE_COLOR_CODE(anRed[iColor], anGreen[iColor],
571 : anBlue[iColor]));
572 0 : panPixelCountPerColorTableEntry->push_back(*pnColor);
573 : }
574 : else
575 : {
576 9 : T *pnColor = HISTOGRAM(histogram, nCLevels, anRed[iColor],
577 : anGreen[iColor], anBlue[iColor]);
578 9 : panPixelCountPerColorTableEntry->push_back(*pnColor);
579 : }
580 : }
581 : }
582 2 : goto end_and_cleanup;
583 : }
584 :
585 : /* ==================================================================== */
586 : /* STEP 3: continually subdivide boxes until no more free */
587 : /* boxes remain or until all colors assigned. */
588 : /* ==================================================================== */
589 7413 : while (freeboxes != nullptr)
590 : {
591 7365 : auto ptr = largest_box(usedboxes);
592 7365 : if (ptr != nullptr)
593 7347 : splitbox(ptr, histogram, psHashHistogram, nCLevels, &freeboxes,
594 : &usedboxes, pabyRedBand, pabyGreenBand, pabyBlueBand,
595 : nPixels);
596 : else
597 18 : freeboxes = nullptr;
598 : }
599 :
600 : /* ==================================================================== */
601 : /* STEP 4: assign colors to all boxes */
602 : /* ==================================================================== */
603 : {
604 48 : Colorbox *ptr = usedboxes;
605 48 : if (panPixelCountPerColorTableEntry)
606 : {
607 20 : panPixelCountPerColorTableEntry->clear();
608 20 : panPixelCountPerColorTableEntry->reserve(nColors);
609 : }
610 7443 : for (int i = 0; ptr != nullptr; ++i, ptr = ptr->next)
611 : {
612 7395 : const GDALColorEntry sEntry = {
613 7395 : static_cast<GByte>(((ptr->rmin + ptr->rmax) << nColorShift) /
614 : 2),
615 7395 : static_cast<GByte>(((ptr->gmin + ptr->gmax) << nColorShift) /
616 : 2),
617 7395 : static_cast<GByte>(((ptr->bmin + ptr->bmax) << nColorShift) /
618 : 2),
619 : 255};
620 7395 : GDALSetColorEntry(hColorTable, i, &sEntry);
621 7395 : if (panPixelCountPerColorTableEntry)
622 798 : panPixelCountPerColorTableEntry->push_back(
623 798 : static_cast<T>(ptr->total));
624 : }
625 : }
626 :
627 48 : end_and_cleanup:
628 50 : CPLFree(pabyRedLine);
629 50 : CPLFree(pabyGreenLine);
630 50 : CPLFree(pabyBlueLine);
631 :
632 : // We're done with the boxes now.
633 50 : CPLFree(box_list);
634 50 : freeboxes = nullptr;
635 50 : usedboxes = nullptr;
636 :
637 50 : if (panHistogram == nullptr)
638 40 : CPLFree(histogram);
639 :
640 50 : return err;
641 : }
642 :
643 : /************************************************************************/
644 : /* largest_box() */
645 : /************************************************************************/
646 :
647 7365 : static Colorbox *largest_box(Colorbox *usedboxes)
648 : {
649 7365 : Colorbox *b = nullptr;
650 :
651 900336 : for (Colorbox *p = usedboxes; p != nullptr; p = p->next)
652 : {
653 892971 : if ((p->rmax > p->rmin || p->gmax > p->gmin || p->bmax > p->bmin) &&
654 732940 : (b == nullptr || p->total > b->total))
655 : {
656 66668 : b = p;
657 : }
658 : }
659 7365 : return b;
660 : }
661 :
662 743 : static void shrinkboxFromBand(Colorbox *ptr, const GByte *pabyRedBand,
663 : const GByte *pabyGreenBand,
664 : const GByte *pabyBlueBand, GUIntBig nPixels)
665 : {
666 743 : int rmin_new = 255;
667 743 : int rmax_new = 0;
668 743 : int gmin_new = 255;
669 743 : int gmax_new = 0;
670 743 : int bmin_new = 255;
671 743 : int bmax_new = 0;
672 35470300 : for (GUIntBig i = 0; i < nPixels; i++)
673 : {
674 35469500 : const int iR = pabyRedBand[i];
675 35469500 : const int iG = pabyGreenBand[i];
676 35469500 : const int iB = pabyBlueBand[i];
677 35469500 : if (iR >= ptr->rmin && iR <= ptr->rmax && iG >= ptr->gmin &&
678 5781920 : iG <= ptr->gmax && iB >= ptr->bmin && iB <= ptr->bmax)
679 : {
680 2154140 : if (iR < rmin_new)
681 3143 : rmin_new = iR;
682 2154140 : if (iR > rmax_new)
683 4274 : rmax_new = iR;
684 2154140 : if (iG < gmin_new)
685 3392 : gmin_new = iG;
686 2154140 : if (iG > gmax_new)
687 4635 : gmax_new = iG;
688 2154140 : if (iB < bmin_new)
689 4455 : bmin_new = iB;
690 2154140 : if (iB > bmax_new)
691 3314 : bmax_new = iB;
692 : }
693 : }
694 :
695 743 : CPLAssert(rmin_new >= ptr->rmin && rmin_new <= rmax_new &&
696 : rmax_new <= ptr->rmax);
697 743 : CPLAssert(gmin_new >= ptr->gmin && gmin_new <= gmax_new &&
698 : gmax_new <= ptr->gmax);
699 743 : CPLAssert(bmin_new >= ptr->bmin && bmin_new <= bmax_new &&
700 : bmax_new <= ptr->bmax);
701 743 : ptr->rmin = rmin_new;
702 743 : ptr->rmax = rmax_new;
703 743 : ptr->gmin = gmin_new;
704 743 : ptr->gmax = gmax_new;
705 743 : ptr->bmin = bmin_new;
706 743 : ptr->bmax = bmax_new;
707 743 : }
708 :
709 3415 : static void shrinkboxFromHashHistogram(Colorbox *box,
710 : const HashHistogram *psHashHistogram)
711 : {
712 3415 : if (box->rmax > box->rmin)
713 : {
714 4917 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
715 : {
716 54728 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
717 : {
718 770529 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
719 : {
720 720718 : if (FindColorCount(psHashHistogram,
721 1441440 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
722 : {
723 3313 : box->rmin = ir;
724 3313 : goto have_rmin;
725 : }
726 : }
727 : }
728 : }
729 : }
730 102 : have_rmin:
731 3415 : if (box->rmax > box->rmin)
732 : {
733 6050 : for (int ir = box->rmax; ir >= box->rmin; --ir)
734 : {
735 65302 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
736 : {
737 1301360 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
738 : {
739 1242110 : if (FindColorCount(psHashHistogram,
740 2484210 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
741 : {
742 3310 : box->rmax = ir;
743 3310 : goto have_rmax;
744 : }
745 : }
746 : }
747 : }
748 : }
749 :
750 105 : have_rmax:
751 3415 : if (box->gmax > box->gmin)
752 : {
753 7027 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
754 : {
755 71453 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
756 : {
757 1233930 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
758 : {
759 1169500 : if (FindColorCount(psHashHistogram,
760 2339010 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
761 : {
762 3252 : box->gmin = ig;
763 3252 : goto have_gmin;
764 : }
765 : }
766 : }
767 : }
768 : }
769 :
770 163 : have_gmin:
771 3415 : if (box->gmax > box->gmin)
772 : {
773 8709 : for (int ig = box->gmax; ig >= box->gmin; --ig)
774 : {
775 99428 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
776 : {
777 93961 : int ib = box->bmin;
778 1785560 : for (; ib <= box->bmax; ++ib)
779 : {
780 1694840 : if (FindColorCount(psHashHistogram,
781 3389690 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
782 : {
783 3242 : box->gmax = ig;
784 3242 : goto have_gmax;
785 : }
786 : }
787 : }
788 : }
789 : }
790 :
791 173 : have_gmax:
792 3415 : if (box->bmax > box->bmin)
793 : {
794 5004 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
795 : {
796 36216 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
797 : {
798 572282 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
799 : {
800 541070 : if (FindColorCount(psHashHistogram,
801 1082140 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
802 : {
803 3312 : box->bmin = ib;
804 3312 : goto have_bmin;
805 : }
806 : }
807 : }
808 : }
809 : }
810 :
811 103 : have_bmin:
812 3415 : if (box->bmax > box->bmin)
813 : {
814 6196 : for (int ib = box->bmax; ib >= box->bmin; --ib)
815 : {
816 48883 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
817 : {
818 675385 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
819 : {
820 632698 : if (FindColorCount(psHashHistogram,
821 1265400 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
822 : {
823 3298 : box->bmax = ib;
824 3298 : goto have_bmax;
825 : }
826 : }
827 : }
828 : }
829 : }
830 :
831 117 : have_bmax:;
832 3415 : }
833 :
834 : /************************************************************************/
835 : /* splitbox() */
836 : /************************************************************************/
837 : template <class T>
838 7347 : static void splitbox(Colorbox *ptr, const T *histogram,
839 : const HashHistogram *psHashHistogram, int nCLevels,
840 : Colorbox **pfreeboxes, Colorbox **pusedboxes,
841 : GByte *pabyRedBand, GByte *pabyGreenBand,
842 : GByte *pabyBlueBand, T nPixels)
843 : {
844 7347 : T hist2[256] = {};
845 7347 : int first = 0;
846 7347 : int last = 0;
847 :
848 : enum
849 : {
850 : RED,
851 : GREEN,
852 : BLUE
853 : } axis;
854 :
855 : // See which axis is the largest, do a histogram along that axis. Split at
856 : // median point. Contract both new boxes to fit points and return.
857 : {
858 7347 : int i = ptr->rmax - ptr->rmin;
859 7347 : if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
860 2466 : axis = RED;
861 4881 : else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
862 2385 : axis = GREEN;
863 : else
864 2496 : axis = BLUE;
865 : }
866 : // Get histogram along longest axis.
867 7347 : const GUInt32 nIters = (ptr->rmax - ptr->rmin + 1) *
868 7347 : (ptr->gmax - ptr->gmin + 1) *
869 7347 : (ptr->bmax - ptr->bmin + 1);
870 :
871 7347 : switch (axis)
872 : {
873 2466 : case RED:
874 : {
875 2466 : if (nPixels != 0 && nIters > nPixels)
876 : {
877 229 : const int rmin = ptr->rmin;
878 229 : const int rmax = ptr->rmax;
879 229 : const int gmin = ptr->gmin;
880 229 : const int gmax = ptr->gmax;
881 229 : const int bmin = ptr->bmin;
882 229 : const int bmax = ptr->bmax;
883 11014200 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
884 : {
885 11014000 : int iR = pabyRedBand[iPixel];
886 11014000 : int iG = pabyGreenBand[iPixel];
887 11014000 : int iB = pabyBlueBand[iPixel];
888 11014000 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
889 1379080 : iB >= bmin && iB <= bmax)
890 : {
891 893690 : hist2[iR]++;
892 : }
893 229 : }
894 : }
895 2237 : else if (psHashHistogram)
896 : {
897 472 : T *histp = &hist2[ptr->rmin];
898 11243 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
899 : {
900 10771 : *histp = 0;
901 226089 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
902 : {
903 3704290 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
904 : {
905 3488980 : *histp += FindColorCount(
906 3488980 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
907 : }
908 : }
909 10771 : histp++;
910 : }
911 : }
912 : else
913 : {
914 1765 : T *histp = &hist2[ptr->rmin];
915 16617 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
916 : {
917 14852 : *histp = 0;
918 188612 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
919 : {
920 : const T *iptr =
921 173760 : HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
922 3151039 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
923 2977279 : *histp += *iptr++;
924 : }
925 14852 : histp++;
926 : }
927 : }
928 2466 : first = ptr->rmin;
929 2466 : last = ptr->rmax;
930 2466 : break;
931 : }
932 2385 : case GREEN:
933 : {
934 2385 : if (nPixels != 0 && nIters > nPixels)
935 : {
936 153 : const int rmin = ptr->rmin;
937 153 : const int rmax = ptr->rmax;
938 153 : const int gmin = ptr->gmin;
939 153 : const int gmax = ptr->gmax;
940 153 : const int bmin = ptr->bmin;
941 153 : const int bmax = ptr->bmax;
942 7442490 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
943 : {
944 7442340 : const int iR = pabyRedBand[iPixel];
945 7442340 : const int iG = pabyGreenBand[iPixel];
946 7442340 : const int iB = pabyBlueBand[iPixel];
947 7442340 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
948 1575360 : iB >= bmin && iB <= bmax)
949 : {
950 948439 : hist2[iG]++;
951 : }
952 153 : }
953 : }
954 2232 : else if (psHashHistogram)
955 : {
956 476 : T *histp = &hist2[ptr->gmin];
957 12600 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
958 : {
959 12124 : *histp = 0;
960 183557 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
961 : {
962 3551710 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
963 : {
964 3380280 : *histp += FindColorCount(
965 3380280 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
966 : }
967 : }
968 12124 : histp++;
969 : }
970 : }
971 : else
972 : {
973 1756 : T *histp = &hist2[ptr->gmin];
974 13924 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
975 : {
976 12168 : *histp = 0;
977 110326 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
978 : {
979 : const T *iptr =
980 98158 : HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
981 1913429 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
982 1815270 : *histp += *iptr++;
983 : }
984 12168 : histp++;
985 : }
986 : }
987 2385 : first = ptr->gmin;
988 2385 : last = ptr->gmax;
989 2385 : break;
990 : }
991 2496 : case BLUE:
992 : {
993 2496 : if (nPixels != 0 && nIters > nPixels)
994 : {
995 199 : const int rmin = ptr->rmin;
996 199 : const int rmax = ptr->rmax;
997 199 : const int gmin = ptr->gmin;
998 199 : const int gmax = ptr->gmax;
999 199 : const int bmin = ptr->bmin;
1000 199 : const int bmax = ptr->bmax;
1001 9614150 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
1002 : {
1003 9613950 : const int iR = pabyRedBand[iPixel];
1004 9613950 : const int iG = pabyGreenBand[iPixel];
1005 9613950 : const int iB = pabyBlueBand[iPixel];
1006 9613950 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
1007 832271 : iB >= bmin && iB <= bmax)
1008 : {
1009 649243 : hist2[iB]++;
1010 : }
1011 199 : }
1012 : }
1013 2297 : else if (psHashHistogram)
1014 : {
1015 573 : T *histp = &hist2[ptr->bmin];
1016 16540 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
1017 : {
1018 15967 : *histp = 0;
1019 252363 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
1020 : {
1021 4915820 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
1022 : {
1023 4679430 : *histp += FindColorCount(
1024 4679430 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
1025 : }
1026 : }
1027 15967 : histp++;
1028 : }
1029 : }
1030 : else
1031 : {
1032 1724 : T *histp = &hist2[ptr->bmin];
1033 15464 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
1034 : {
1035 13740 : *histp = 0;
1036 136611 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
1037 : {
1038 : const T *iptr =
1039 122871 : HISTOGRAM(histogram, nCLevels, ir, ptr->gmin, ib);
1040 2358380 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
1041 : {
1042 2235514 : *histp += *iptr;
1043 2235514 : iptr += nCLevels;
1044 : }
1045 : }
1046 13740 : histp++;
1047 : }
1048 : }
1049 2496 : first = ptr->bmin;
1050 2496 : last = ptr->bmax;
1051 2496 : break;
1052 : }
1053 : }
1054 : // Find median point.
1055 7347 : T *histp = &hist2[first];
1056 7347 : int i = first; // TODO(schwehr): Rename i.
1057 : {
1058 7347 : T sum = 0;
1059 7347 : T sum2 = static_cast<T>(ptr->total / 2);
1060 57172 : for (; i <= last && (sum += *histp++) < sum2; ++i)
1061 : {
1062 : }
1063 : }
1064 7347 : if (i == first)
1065 1799 : i++;
1066 :
1067 : // Create new box, re-allocate points.
1068 7347 : Colorbox *new_cb = *pfreeboxes;
1069 7347 : *pfreeboxes = new_cb->next;
1070 7347 : if (*pfreeboxes)
1071 7317 : (*pfreeboxes)->prev = nullptr;
1072 7347 : if (*pusedboxes)
1073 7347 : (*pusedboxes)->prev = new_cb;
1074 7347 : new_cb->next = *pusedboxes;
1075 7347 : *pusedboxes = new_cb;
1076 :
1077 7347 : histp = &hist2[first];
1078 : {
1079 7347 : T sum1 = 0;
1080 58971 : for (int j = first; j < i; j++)
1081 51624 : sum1 += *histp++;
1082 7347 : T sum2 = 0;
1083 91652 : for (int j = i; j <= last; j++)
1084 84305 : sum2 += *histp++;
1085 7347 : new_cb->total = sum1;
1086 7347 : ptr->total = sum2;
1087 : }
1088 7347 : new_cb->rmin = ptr->rmin;
1089 7347 : new_cb->rmax = ptr->rmax;
1090 7347 : new_cb->gmin = ptr->gmin;
1091 7347 : new_cb->gmax = ptr->gmax;
1092 7347 : new_cb->bmin = ptr->bmin;
1093 7347 : new_cb->bmax = ptr->bmax;
1094 7347 : switch (axis)
1095 : {
1096 2466 : case RED:
1097 2466 : new_cb->rmax = i - 1;
1098 2466 : ptr->rmin = i;
1099 2466 : break;
1100 2385 : case GREEN:
1101 2385 : new_cb->gmax = i - 1;
1102 2385 : ptr->gmin = i;
1103 2385 : break;
1104 2496 : case BLUE:
1105 2496 : new_cb->bmax = i - 1;
1106 2496 : ptr->bmin = i;
1107 2496 : break;
1108 : }
1109 :
1110 7347 : if (nPixels != 0 &&
1111 2295 : static_cast<T>(new_cb->rmax - new_cb->rmin + 1) *
1112 2295 : static_cast<T>(new_cb->gmax - new_cb->gmin + 1) *
1113 2295 : static_cast<T>(new_cb->bmax - new_cb->bmin + 1) >
1114 : nPixels)
1115 : {
1116 325 : shrinkboxFromBand(new_cb, pabyRedBand, pabyGreenBand, pabyBlueBand,
1117 : nPixels);
1118 : }
1119 7022 : else if (psHashHistogram != nullptr)
1120 : {
1121 1748 : shrinkboxFromHashHistogram(new_cb, psHashHistogram);
1122 : }
1123 : else
1124 : {
1125 5274 : shrinkbox(new_cb, histogram, nCLevels);
1126 : }
1127 :
1128 7347 : if (nPixels != 0 && static_cast<T>(ptr->rmax - ptr->rmin + 1) *
1129 2295 : static_cast<T>(ptr->gmax - ptr->gmin + 1) *
1130 2295 : static_cast<T>(ptr->bmax - ptr->bmin + 1) >
1131 : nPixels)
1132 : {
1133 418 : shrinkboxFromBand(ptr, pabyRedBand, pabyGreenBand, pabyBlueBand,
1134 : nPixels);
1135 : }
1136 6929 : else if (psHashHistogram != nullptr)
1137 : {
1138 1667 : shrinkboxFromHashHistogram(ptr, psHashHistogram);
1139 : }
1140 : else
1141 : {
1142 5262 : shrinkbox(ptr, histogram, nCLevels);
1143 : }
1144 7347 : }
1145 :
1146 : /************************************************************************/
1147 : /* shrinkbox() */
1148 : /************************************************************************/
1149 : template <class T>
1150 10536 : static void shrinkbox(Colorbox *box, const T *histogram, int nCLevels)
1151 : {
1152 10536 : if (box->rmax > box->rmin)
1153 : {
1154 7042 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1155 : {
1156 20590 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1157 : {
1158 : const T *histp =
1159 19589 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1160 244021 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1161 : {
1162 230473 : if (*histp++ != 0)
1163 : {
1164 6041 : box->rmin = ir;
1165 6041 : goto have_rmin;
1166 : }
1167 : }
1168 : }
1169 : }
1170 : }
1171 4495 : have_rmin:
1172 10536 : if (box->rmax > box->rmin)
1173 : {
1174 7820 : for (int ir = box->rmax; ir >= box->rmin; --ir)
1175 : {
1176 34920 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1177 : {
1178 : const T *histp =
1179 33087 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1180 463269 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1181 : {
1182 436169 : if (*histp++ != 0)
1183 : {
1184 5987 : box->rmax = ir;
1185 5987 : goto have_rmax;
1186 : }
1187 : }
1188 : }
1189 : }
1190 : }
1191 :
1192 4549 : have_rmax:
1193 10536 : if (box->gmax > box->gmin)
1194 : {
1195 9144 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1196 : {
1197 29845 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1198 : {
1199 : const T *histp =
1200 27743 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1201 401339 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1202 : {
1203 380620 : if (*histp++ != 0)
1204 : {
1205 7024 : box->gmin = ig;
1206 7024 : goto have_gmin;
1207 : }
1208 : }
1209 : }
1210 : }
1211 : }
1212 :
1213 3494 : have_gmin:
1214 10536 : if (box->gmax > box->gmin)
1215 : {
1216 11393 : for (int ig = box->gmax; ig >= box->gmin; --ig)
1217 : {
1218 51896 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1219 : {
1220 : const T *histp =
1221 47405 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1222 666264 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1223 : {
1224 625743 : if (*histp++ != 0)
1225 : {
1226 6884 : box->gmax = ig;
1227 6884 : goto have_gmax;
1228 : }
1229 : }
1230 : }
1231 : }
1232 : }
1233 :
1234 3634 : have_gmax:
1235 10536 : if (box->bmax > box->bmin)
1236 : {
1237 7845 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1238 : {
1239 18296 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1240 : {
1241 : const T *histp =
1242 17082 : HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
1243 166256 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1244 : {
1245 155789 : if (*histp != 0)
1246 : {
1247 6615 : box->bmin = ib;
1248 6615 : goto have_bmin;
1249 : }
1250 149174 : histp += nCLevels;
1251 : }
1252 : }
1253 : }
1254 : }
1255 :
1256 3905 : have_bmin:
1257 10536 : if (box->bmax > box->bmin)
1258 : {
1259 10612 : for (int ib = box->bmax; ib >= box->bmin; --ib)
1260 : {
1261 40412 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1262 : {
1263 : const T *histp =
1264 36299 : HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
1265 379858 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1266 : {
1267 350042 : if (*histp != 0)
1268 : {
1269 6483 : box->bmax = ib;
1270 6483 : goto have_bmax;
1271 : }
1272 343559 : histp += nCLevels;
1273 : }
1274 : }
1275 : }
1276 : }
1277 :
1278 4037 : have_bmax:;
1279 10536 : }
1280 :
1281 : // Explicitly instantiate template functions.
1282 : template int GDALComputeMedianCutPCTInternal<GUInt32>(
1283 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1284 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1285 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1286 : GUInt32 *panHistogram, GDALColorTableH hColorTable,
1287 : GDALProgressFunc pfnProgress, void *pProgressArg,
1288 : std::vector<GUInt32> *panPixelCountPerColorTableEntry);
1289 :
1290 : template int GDALComputeMedianCutPCTInternal<GUIntBig>(
1291 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1292 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1293 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1294 : GUIntBig *panHistogram, GDALColorTableH hColorTable,
1295 : GDALProgressFunc pfnProgress, void *pProgressArg,
1296 : std::vector<GUIntBig> *panPixelCountPerColorTableEntry);
1297 :
1298 0 : int GDALComputeMedianCutPCT(
1299 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1300 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1301 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1302 : GUInt32 *panHistogram, GDALColorTableH hColorTable,
1303 : GDALProgressFunc pfnProgress, void *pProgressArg,
1304 : std::vector<GUInt32> *panPixelCountPerColorTableEntry)
1305 : {
1306 0 : return GDALComputeMedianCutPCTInternal<GUInt32>(
1307 : hRed, hGreen, hBlue, pabyRedBand, pabyGreenBand, pabyBlueBand,
1308 : pfnIncludePixel, nColors, nBits, panHistogram, hColorTable, pfnProgress,
1309 0 : pProgressArg, panPixelCountPerColorTableEntry);
1310 : }
1311 :
1312 21 : int GDALComputeMedianCutPCT(
1313 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1314 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1315 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1316 : GUIntBig *panHistogram, GDALColorTableH hColorTable,
1317 : GDALProgressFunc pfnProgress, void *pProgressArg,
1318 : std::vector<GUIntBig> *panPixelCountPerColorTableEntry)
1319 : {
1320 21 : return GDALComputeMedianCutPCTInternal<GUIntBig>(
1321 : hRed, hGreen, hBlue, pabyRedBand, pabyGreenBand, pabyBlueBand,
1322 : pfnIncludePixel, nColors, nBits, panHistogram, hColorTable, pfnProgress,
1323 21 : pProgressArg, panPixelCountPerColorTableEntry);
1324 : }
|
