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 1391663 : template <typename T> static T *HISTOGRAM(T *h, int n, int r, int g, int b)
42 : {
43 1391663 : const int index = (r * n + g) * n + b;
44 1391663 : 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 27 : 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)
276 :
277 : {
278 27 : VALIDATE_POINTER1(hRed, "GDALComputeMedianCutPCT", CE_Failure);
279 27 : VALIDATE_POINTER1(hGreen, "GDALComputeMedianCutPCT", CE_Failure);
280 27 : VALIDATE_POINTER1(hBlue, "GDALComputeMedianCutPCT", CE_Failure);
281 :
282 27 : CPLErr err = CE_None;
283 :
284 : /* -------------------------------------------------------------------- */
285 : /* Validate parameters. */
286 : /* -------------------------------------------------------------------- */
287 27 : const int nXSize = GDALGetRasterBandXSize(hRed);
288 27 : const int nYSize = GDALGetRasterBandYSize(hRed);
289 :
290 27 : if (GDALGetRasterBandXSize(hGreen) != nXSize ||
291 27 : GDALGetRasterBandYSize(hGreen) != nYSize ||
292 81 : GDALGetRasterBandXSize(hBlue) != nXSize ||
293 27 : 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 27 : 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 27 : 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 27 : 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 27 : if (pfnProgress == nullptr)
329 21 : pfnProgress = GDALDummyProgress;
330 :
331 27 : int nSrcNoData = -1;
332 27 : GDALRasterBandH hMaskBand = nullptr;
333 27 : int bSrcHasNoDataR = FALSE;
334 27 : const double dfSrcNoDataR = GDALGetRasterNoDataValue(hRed, &bSrcHasNoDataR);
335 27 : int bSrcHasNoDataG = FALSE;
336 : const double dfSrcNoDataG =
337 27 : GDALGetRasterNoDataValue(hGreen, &bSrcHasNoDataG);
338 27 : int bSrcHasNoDataB = FALSE;
339 : const double dfSrcNoDataB =
340 27 : GDALGetRasterNoDataValue(hBlue, &bSrcHasNoDataB);
341 27 : 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 25 : const int nMaskFlags = GDALGetMaskFlags(hRed);
351 25 : if ((nMaskFlags & GMF_PER_DATASET))
352 1 : hMaskBand = GDALGetMaskBand(hRed);
353 : }
354 :
355 : /* ==================================================================== */
356 : /* STEP 1: create empty boxes. */
357 : /* ==================================================================== */
358 27 : if (static_cast<GUInt32>(nXSize) >
359 27 : 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 27 : T nPixels = 0;
367 10 : if (nBits == 8 && pabyRedBand != nullptr && pabyGreenBand != nullptr &&
368 37 : 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 27 : const int nCLevels = 1 << nBits;
376 27 : const int nCLevelsCube = nCLevels * nCLevels * nCLevels;
377 27 : T *histogram = nullptr;
378 27 : HashHistogram *psHashHistogram = nullptr;
379 27 : 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 17 : static_cast<T *>(VSI_CALLOC_VERBOSE(nCLevelsCube, sizeof(T)));
401 17 : if (histogram == nullptr)
402 : {
403 0 : return CE_Failure;
404 : }
405 : }
406 : Colorbox *box_list =
407 27 : static_cast<Colorbox *>(CPLMalloc(nColors * sizeof(Colorbox)));
408 27 : Colorbox *freeboxes = box_list;
409 27 : freeboxes[0].next = &freeboxes[1];
410 27 : freeboxes[0].prev = nullptr;
411 6394 : for (int i = 1; i < nColors - 1; ++i)
412 : {
413 6367 : freeboxes[i].next = &freeboxes[i + 1];
414 6367 : freeboxes[i].prev = &freeboxes[i - 1];
415 : }
416 27 : freeboxes[nColors - 1].next = nullptr;
417 27 : freeboxes[nColors - 1].prev = &freeboxes[nColors - 2];
418 :
419 : /* ==================================================================== */
420 : /* Build histogram. */
421 : /* ==================================================================== */
422 :
423 : /* -------------------------------------------------------------------- */
424 : /* Initialize the box datastructures. */
425 : /* -------------------------------------------------------------------- */
426 27 : Colorbox *freeboxes_before = freeboxes;
427 27 : freeboxes = freeboxes_before->next;
428 27 : if (freeboxes)
429 27 : freeboxes->prev = nullptr;
430 :
431 27 : Colorbox *usedboxes = freeboxes_before;
432 27 : usedboxes->next = nullptr;
433 27 : usedboxes->rmin = 999;
434 27 : usedboxes->gmin = 999;
435 27 : usedboxes->bmin = 999;
436 27 : usedboxes->rmax = -1;
437 27 : usedboxes->gmax = -1;
438 27 : usedboxes->bmax = -1;
439 27 : usedboxes->total =
440 27 : 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 27 : const int nColorShift = 8 - nBits;
448 27 : int nColorCounter = 0;
449 27 : GByte anRed[256] = {};
450 27 : GByte anGreen[256] = {};
451 27 : GByte anBlue[256] = {};
452 :
453 27 : GByte *pabyRedLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
454 27 : GByte *pabyGreenLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
455 27 : GByte *pabyBlueLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
456 0 : std::unique_ptr<GByte, VSIFreeReleaser> pabyMask;
457 27 : if (hMaskBand)
458 1 : pabyMask.reset(static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)));
459 :
460 27 : if (pabyRedLine == nullptr || pabyGreenLine == nullptr ||
461 54 : pabyBlueLine == nullptr || (hMaskBand && !pabyMask))
462 : {
463 0 : err = CE_Failure;
464 0 : goto end_and_cleanup;
465 : }
466 :
467 4651 : for (int iLine = 0; iLine < nYSize; iLine++)
468 : {
469 4624 : 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 4624 : err = GDALRasterIO(hRed, GF_Read, 0, iLine, nXSize, 1, pabyRedLine,
478 : nXSize, 1, GDT_UInt8, 0, 0);
479 4624 : if (err == CE_None)
480 4624 : err = GDALRasterIO(hGreen, GF_Read, 0, iLine, nXSize, 1,
481 : pabyGreenLine, nXSize, 1, GDT_UInt8, 0, 0);
482 4624 : if (err == CE_None)
483 4624 : err = GDALRasterIO(hBlue, GF_Read, 0, iLine, nXSize, 1,
484 : pabyBlueLine, nXSize, 1, GDT_UInt8, 0, 0);
485 4624 : if (err == CE_None && hMaskBand)
486 150 : err = GDALRasterIO(hMaskBand, GF_Read, 0, iLine, nXSize, 1,
487 150 : pabyMask.get(), nXSize, 1, GDT_UInt8, 0, 0);
488 4624 : if (err != CE_None)
489 0 : goto end_and_cleanup;
490 :
491 1357170 : for (int iPixel = 0; iPixel < nXSize; iPixel++)
492 : {
493 2783730 : if ((pabyRedLine[iPixel] == nSrcNoData &&
494 33554 : pabyGreenLine[iPixel] == nSrcNoData &&
495 2729390 : pabyBlueLine[iPixel] == nSrcNoData) ||
496 1343540 : (pabyMask && (pabyMask.get())[iPixel] == 0))
497 : {
498 45088 : continue;
499 : }
500 :
501 1307460 : const int nRed = pabyRedLine[iPixel] >> nColorShift;
502 1307460 : const int nGreen = pabyGreenLine[iPixel] >> nColorShift;
503 1307460 : const int nBlue = pabyBlueLine[iPixel] >> nColorShift;
504 :
505 1307460 : usedboxes->rmin = std::min(usedboxes->rmin, nRed);
506 1307460 : usedboxes->gmin = std::min(usedboxes->gmin, nGreen);
507 1307460 : usedboxes->bmin = std::min(usedboxes->bmin, nBlue);
508 1307460 : usedboxes->rmax = std::max(usedboxes->rmax, nRed);
509 1307460 : usedboxes->gmax = std::max(usedboxes->gmax, nGreen);
510 1307460 : usedboxes->bmax = std::max(usedboxes->bmax, nBlue);
511 :
512 : bool bFirstOccurrence;
513 1307460 : 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 917812 : HISTOGRAM(histogram, nCLevels, nRed, nGreen, nBlue);
524 917812 : bFirstOccurrence = (*pnColor == 0);
525 917812 : (*pnColor)++;
526 : }
527 1307460 : if (bFirstOccurrence)
528 : {
529 104280 : if (nColorShift == 0 && nColorCounter < nColors)
530 : {
531 2305 : anRed[nColorCounter] = static_cast<GByte>(nRed);
532 2305 : anGreen[nColorCounter] = static_cast<GByte>(nGreen);
533 2305 : anBlue[nColorCounter] = static_cast<GByte>(nBlue);
534 : }
535 104280 : nColorCounter++;
536 : }
537 : }
538 : }
539 :
540 27 : 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 27 : if (nColorShift == 0 && nColorCounter <= nColors)
548 : {
549 : #if DEBUG_VERBOSE
550 : CPLDebug("MEDIAN_CUT", "%d colors found <= %d", nColorCounter, nColors);
551 : #endif
552 2 : for (int iColor = 0; iColor < nColorCounter; iColor++)
553 : {
554 1 : const GDALColorEntry sEntry = {static_cast<GByte>(anRed[iColor]),
555 1 : static_cast<GByte>(anGreen[iColor]),
556 1 : static_cast<GByte>(anBlue[iColor]),
557 : 255};
558 1 : GDALSetColorEntry(hColorTable, iColor, &sEntry);
559 : }
560 1 : goto end_and_cleanup;
561 : }
562 :
563 : /* ==================================================================== */
564 : /* STEP 3: continually subdivide boxes until no more free */
565 : /* boxes remain or until all colors assigned. */
566 : /* ==================================================================== */
567 6165 : while (freeboxes != nullptr)
568 : {
569 6139 : auto ptr = largest_box(usedboxes);
570 6139 : if (ptr != nullptr)
571 6139 : splitbox(ptr, histogram, psHashHistogram, nCLevels, &freeboxes,
572 : &usedboxes, pabyRedBand, pabyGreenBand, pabyBlueBand,
573 : nPixels);
574 : else
575 0 : freeboxes = nullptr;
576 : }
577 :
578 : /* ==================================================================== */
579 : /* STEP 4: assign colors to all boxes */
580 : /* ==================================================================== */
581 : {
582 26 : Colorbox *ptr = usedboxes;
583 6191 : for (int i = 0; ptr != nullptr; ++i, ptr = ptr->next)
584 : {
585 6165 : const GDALColorEntry sEntry = {
586 6165 : static_cast<GByte>(((ptr->rmin + ptr->rmax) << nColorShift) /
587 : 2),
588 6165 : static_cast<GByte>(((ptr->gmin + ptr->gmax) << nColorShift) /
589 : 2),
590 6165 : static_cast<GByte>(((ptr->bmin + ptr->bmax) << nColorShift) /
591 : 2),
592 : 255};
593 6165 : GDALSetColorEntry(hColorTable, i, &sEntry);
594 : }
595 : }
596 :
597 26 : end_and_cleanup:
598 27 : CPLFree(pabyRedLine);
599 27 : CPLFree(pabyGreenLine);
600 27 : CPLFree(pabyBlueLine);
601 :
602 : // We're done with the boxes now.
603 27 : CPLFree(box_list);
604 27 : freeboxes = nullptr;
605 27 : usedboxes = nullptr;
606 :
607 27 : if (panHistogram == nullptr)
608 17 : CPLFree(histogram);
609 :
610 27 : return err;
611 : }
612 :
613 : /************************************************************************/
614 : /* largest_box() */
615 : /************************************************************************/
616 :
617 6139 : static Colorbox *largest_box(Colorbox *usedboxes)
618 : {
619 6139 : Colorbox *b = nullptr;
620 :
621 788882 : for (Colorbox *p = usedboxes; p != nullptr; p = p->next)
622 : {
623 782743 : if ((p->rmax > p->rmin || p->gmax > p->gmin || p->bmax > p->bmin) &&
624 659765 : (b == nullptr || p->total > b->total))
625 : {
626 61013 : b = p;
627 : }
628 : }
629 6139 : return b;
630 : }
631 :
632 743 : static void shrinkboxFromBand(Colorbox *ptr, const GByte *pabyRedBand,
633 : const GByte *pabyGreenBand,
634 : const GByte *pabyBlueBand, GUIntBig nPixels)
635 : {
636 743 : int rmin_new = 255;
637 743 : int rmax_new = 0;
638 743 : int gmin_new = 255;
639 743 : int gmax_new = 0;
640 743 : int bmin_new = 255;
641 743 : int bmax_new = 0;
642 35470300 : for (GUIntBig i = 0; i < nPixels; i++)
643 : {
644 35469500 : const int iR = pabyRedBand[i];
645 35469500 : const int iG = pabyGreenBand[i];
646 35469500 : const int iB = pabyBlueBand[i];
647 35469500 : if (iR >= ptr->rmin && iR <= ptr->rmax && iG >= ptr->gmin &&
648 5781920 : iG <= ptr->gmax && iB >= ptr->bmin && iB <= ptr->bmax)
649 : {
650 2154140 : if (iR < rmin_new)
651 3143 : rmin_new = iR;
652 2154140 : if (iR > rmax_new)
653 4274 : rmax_new = iR;
654 2154140 : if (iG < gmin_new)
655 3392 : gmin_new = iG;
656 2154140 : if (iG > gmax_new)
657 4635 : gmax_new = iG;
658 2154140 : if (iB < bmin_new)
659 4455 : bmin_new = iB;
660 2154140 : if (iB > bmax_new)
661 3314 : bmax_new = iB;
662 : }
663 : }
664 :
665 743 : CPLAssert(rmin_new >= ptr->rmin && rmin_new <= rmax_new &&
666 : rmax_new <= ptr->rmax);
667 743 : CPLAssert(gmin_new >= ptr->gmin && gmin_new <= gmax_new &&
668 : gmax_new <= ptr->gmax);
669 743 : CPLAssert(bmin_new >= ptr->bmin && bmin_new <= bmax_new &&
670 : bmax_new <= ptr->bmax);
671 743 : ptr->rmin = rmin_new;
672 743 : ptr->rmax = rmax_new;
673 743 : ptr->gmin = gmin_new;
674 743 : ptr->gmax = gmax_new;
675 743 : ptr->bmin = bmin_new;
676 743 : ptr->bmax = bmax_new;
677 743 : }
678 :
679 3415 : static void shrinkboxFromHashHistogram(Colorbox *box,
680 : const HashHistogram *psHashHistogram)
681 : {
682 3415 : if (box->rmax > box->rmin)
683 : {
684 4917 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
685 : {
686 54728 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
687 : {
688 770529 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
689 : {
690 720718 : if (FindColorCount(psHashHistogram,
691 1441440 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
692 : {
693 3313 : box->rmin = ir;
694 3313 : goto have_rmin;
695 : }
696 : }
697 : }
698 : }
699 : }
700 102 : have_rmin:
701 3415 : if (box->rmax > box->rmin)
702 : {
703 6050 : for (int ir = box->rmax; ir >= box->rmin; --ir)
704 : {
705 65302 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
706 : {
707 1301360 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
708 : {
709 1242110 : if (FindColorCount(psHashHistogram,
710 2484210 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
711 : {
712 3310 : box->rmax = ir;
713 3310 : goto have_rmax;
714 : }
715 : }
716 : }
717 : }
718 : }
719 :
720 105 : have_rmax:
721 3415 : if (box->gmax > box->gmin)
722 : {
723 7027 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
724 : {
725 71453 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
726 : {
727 1233930 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
728 : {
729 1169500 : if (FindColorCount(psHashHistogram,
730 2339010 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
731 : {
732 3252 : box->gmin = ig;
733 3252 : goto have_gmin;
734 : }
735 : }
736 : }
737 : }
738 : }
739 :
740 163 : have_gmin:
741 3415 : if (box->gmax > box->gmin)
742 : {
743 8709 : for (int ig = box->gmax; ig >= box->gmin; --ig)
744 : {
745 99428 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
746 : {
747 93961 : int ib = box->bmin;
748 1785560 : for (; ib <= box->bmax; ++ib)
749 : {
750 1694840 : if (FindColorCount(psHashHistogram,
751 3389690 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
752 : {
753 3242 : box->gmax = ig;
754 3242 : goto have_gmax;
755 : }
756 : }
757 : }
758 : }
759 : }
760 :
761 173 : have_gmax:
762 3415 : if (box->bmax > box->bmin)
763 : {
764 5004 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
765 : {
766 36216 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
767 : {
768 572282 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
769 : {
770 541070 : if (FindColorCount(psHashHistogram,
771 1082140 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
772 : {
773 3312 : box->bmin = ib;
774 3312 : goto have_bmin;
775 : }
776 : }
777 : }
778 : }
779 : }
780 :
781 103 : have_bmin:
782 3415 : if (box->bmax > box->bmin)
783 : {
784 6196 : for (int ib = box->bmax; ib >= box->bmin; --ib)
785 : {
786 48883 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
787 : {
788 675385 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
789 : {
790 632698 : if (FindColorCount(psHashHistogram,
791 1265400 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
792 : {
793 3298 : box->bmax = ib;
794 3298 : goto have_bmax;
795 : }
796 : }
797 : }
798 : }
799 : }
800 :
801 117 : have_bmax:;
802 3415 : }
803 :
804 : /************************************************************************/
805 : /* splitbox() */
806 : /************************************************************************/
807 : template <class T>
808 6139 : static void splitbox(Colorbox *ptr, const T *histogram,
809 : const HashHistogram *psHashHistogram, int nCLevels,
810 : Colorbox **pfreeboxes, Colorbox **pusedboxes,
811 : GByte *pabyRedBand, GByte *pabyGreenBand,
812 : GByte *pabyBlueBand, T nPixels)
813 : {
814 6139 : T hist2[256] = {};
815 6139 : int first = 0;
816 6139 : int last = 0;
817 :
818 : enum
819 : {
820 : RED,
821 : GREEN,
822 : BLUE
823 : } axis;
824 :
825 : // See which axis is the largest, do a histogram along that axis. Split at
826 : // median point. Contract both new boxes to fit points and return.
827 : {
828 6139 : int i = ptr->rmax - ptr->rmin;
829 6139 : if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
830 2001 : axis = RED;
831 4138 : else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
832 2031 : axis = GREEN;
833 : else
834 2107 : axis = BLUE;
835 : }
836 : // Get histogram along longest axis.
837 6139 : const GUInt32 nIters = (ptr->rmax - ptr->rmin + 1) *
838 6139 : (ptr->gmax - ptr->gmin + 1) *
839 6139 : (ptr->bmax - ptr->bmin + 1);
840 :
841 6139 : switch (axis)
842 : {
843 2001 : case RED:
844 : {
845 2001 : if (nPixels != 0 && nIters > nPixels)
846 : {
847 229 : const int rmin = ptr->rmin;
848 229 : const int rmax = ptr->rmax;
849 229 : const int gmin = ptr->gmin;
850 229 : const int gmax = ptr->gmax;
851 229 : const int bmin = ptr->bmin;
852 229 : const int bmax = ptr->bmax;
853 11014200 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
854 : {
855 11014000 : int iR = pabyRedBand[iPixel];
856 11014000 : int iG = pabyGreenBand[iPixel];
857 11014000 : int iB = pabyBlueBand[iPixel];
858 11014000 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
859 1379080 : iB >= bmin && iB <= bmax)
860 : {
861 893690 : hist2[iR]++;
862 : }
863 229 : }
864 : }
865 1772 : else if (psHashHistogram)
866 : {
867 472 : T *histp = &hist2[ptr->rmin];
868 11243 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
869 : {
870 10771 : *histp = 0;
871 226089 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
872 : {
873 3704290 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
874 : {
875 3488980 : *histp += FindColorCount(
876 3488980 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
877 : }
878 : }
879 10771 : histp++;
880 : }
881 : }
882 : else
883 : {
884 1300 : T *histp = &hist2[ptr->rmin];
885 12581 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
886 : {
887 11281 : *histp = 0;
888 144873 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
889 : {
890 : const T *iptr =
891 133592 : HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
892 2460750 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
893 2327160 : *histp += *iptr++;
894 : }
895 11281 : histp++;
896 : }
897 : }
898 2001 : first = ptr->rmin;
899 2001 : last = ptr->rmax;
900 2001 : break;
901 : }
902 2031 : case GREEN:
903 : {
904 2031 : if (nPixels != 0 && nIters > nPixels)
905 : {
906 153 : const int rmin = ptr->rmin;
907 153 : const int rmax = ptr->rmax;
908 153 : const int gmin = ptr->gmin;
909 153 : const int gmax = ptr->gmax;
910 153 : const int bmin = ptr->bmin;
911 153 : const int bmax = ptr->bmax;
912 7442490 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
913 : {
914 7442340 : const int iR = pabyRedBand[iPixel];
915 7442340 : const int iG = pabyGreenBand[iPixel];
916 7442340 : const int iB = pabyBlueBand[iPixel];
917 7442340 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
918 1575360 : iB >= bmin && iB <= bmax)
919 : {
920 948439 : hist2[iG]++;
921 : }
922 153 : }
923 : }
924 1878 : else if (psHashHistogram)
925 : {
926 476 : T *histp = &hist2[ptr->gmin];
927 12600 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
928 : {
929 12124 : *histp = 0;
930 183557 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
931 : {
932 3551710 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
933 : {
934 3380280 : *histp += FindColorCount(
935 3380280 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
936 : }
937 : }
938 12124 : histp++;
939 : }
940 : }
941 : else
942 : {
943 1402 : T *histp = &hist2[ptr->gmin];
944 11707 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
945 : {
946 10305 : *histp = 0;
947 97506 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
948 : {
949 : const T *iptr =
950 87201 : HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
951 1695390 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
952 1608190 : *histp += *iptr++;
953 : }
954 10305 : histp++;
955 : }
956 : }
957 2031 : first = ptr->gmin;
958 2031 : last = ptr->gmax;
959 2031 : break;
960 : }
961 2107 : case BLUE:
962 : {
963 2107 : if (nPixels != 0 && nIters > nPixels)
964 : {
965 199 : const int rmin = ptr->rmin;
966 199 : const int rmax = ptr->rmax;
967 199 : const int gmin = ptr->gmin;
968 199 : const int gmax = ptr->gmax;
969 199 : const int bmin = ptr->bmin;
970 199 : const int bmax = ptr->bmax;
971 9614150 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
972 : {
973 9613950 : const int iR = pabyRedBand[iPixel];
974 9613950 : const int iG = pabyGreenBand[iPixel];
975 9613950 : const int iB = pabyBlueBand[iPixel];
976 9613950 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
977 832271 : iB >= bmin && iB <= bmax)
978 : {
979 649243 : hist2[iB]++;
980 : }
981 199 : }
982 : }
983 1908 : else if (psHashHistogram)
984 : {
985 573 : T *histp = &hist2[ptr->bmin];
986 16540 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
987 : {
988 15967 : *histp = 0;
989 252363 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
990 : {
991 4915820 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
992 : {
993 4679430 : *histp += FindColorCount(
994 4679430 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
995 : }
996 : }
997 15967 : histp++;
998 : }
999 : }
1000 : else
1001 : {
1002 1335 : T *histp = &hist2[ptr->bmin];
1003 12782 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
1004 : {
1005 11447 : *histp = 0;
1006 120522 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
1007 : {
1008 : const T *iptr =
1009 109075 : HISTOGRAM(histogram, nCLevels, ir, ptr->gmin, ib);
1010 2167660 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
1011 : {
1012 2058580 : *histp += *iptr;
1013 2058580 : iptr += nCLevels;
1014 : }
1015 : }
1016 11447 : histp++;
1017 : }
1018 : }
1019 2107 : first = ptr->bmin;
1020 2107 : last = ptr->bmax;
1021 2107 : break;
1022 : }
1023 : }
1024 : // Find median point.
1025 6139 : T *histp = &hist2[first];
1026 6139 : int i = first; // TODO(schwehr): Rename i.
1027 : {
1028 6139 : T sum = 0;
1029 6139 : T sum2 = static_cast<T>(ptr->total / 2);
1030 53761 : for (; i <= last && (sum += *histp++) < sum2; ++i)
1031 : {
1032 : }
1033 : }
1034 6139 : if (i == first)
1035 1347 : i++;
1036 :
1037 : // Create new box, re-allocate points.
1038 6139 : Colorbox *new_cb = *pfreeboxes;
1039 6139 : *pfreeboxes = new_cb->next;
1040 6139 : if (*pfreeboxes)
1041 6113 : (*pfreeboxes)->prev = nullptr;
1042 6139 : if (*pusedboxes)
1043 6139 : (*pusedboxes)->prev = new_cb;
1044 6139 : new_cb->next = *pusedboxes;
1045 6139 : *pusedboxes = new_cb;
1046 :
1047 6139 : histp = &hist2[first];
1048 : {
1049 6139 : T sum1 = 0;
1050 55108 : for (int j = first; j < i; j++)
1051 48969 : sum1 += *histp++;
1052 6139 : T sum2 = 0;
1053 85372 : for (int j = i; j <= last; j++)
1054 79233 : sum2 += *histp++;
1055 6139 : new_cb->total = sum1;
1056 6139 : ptr->total = sum2;
1057 : }
1058 6139 : new_cb->rmin = ptr->rmin;
1059 6139 : new_cb->rmax = ptr->rmax;
1060 6139 : new_cb->gmin = ptr->gmin;
1061 6139 : new_cb->gmax = ptr->gmax;
1062 6139 : new_cb->bmin = ptr->bmin;
1063 6139 : new_cb->bmax = ptr->bmax;
1064 6139 : switch (axis)
1065 : {
1066 2001 : case RED:
1067 2001 : new_cb->rmax = i - 1;
1068 2001 : ptr->rmin = i;
1069 2001 : break;
1070 2031 : case GREEN:
1071 2031 : new_cb->gmax = i - 1;
1072 2031 : ptr->gmin = i;
1073 2031 : break;
1074 2107 : case BLUE:
1075 2107 : new_cb->bmax = i - 1;
1076 2107 : ptr->bmin = i;
1077 2107 : break;
1078 : }
1079 :
1080 6139 : if (nPixels != 0 &&
1081 2295 : static_cast<T>(new_cb->rmax - new_cb->rmin + 1) *
1082 2295 : static_cast<T>(new_cb->gmax - new_cb->gmin + 1) *
1083 2295 : static_cast<T>(new_cb->bmax - new_cb->bmin + 1) >
1084 : nPixels)
1085 : {
1086 325 : shrinkboxFromBand(new_cb, pabyRedBand, pabyGreenBand, pabyBlueBand,
1087 : nPixels);
1088 : }
1089 5814 : else if (psHashHistogram != nullptr)
1090 : {
1091 1748 : shrinkboxFromHashHistogram(new_cb, psHashHistogram);
1092 : }
1093 : else
1094 : {
1095 4066 : shrinkbox(new_cb, histogram, nCLevels);
1096 : }
1097 :
1098 6139 : if (nPixels != 0 && static_cast<T>(ptr->rmax - ptr->rmin + 1) *
1099 2295 : static_cast<T>(ptr->gmax - ptr->gmin + 1) *
1100 2295 : static_cast<T>(ptr->bmax - ptr->bmin + 1) >
1101 : nPixels)
1102 : {
1103 418 : shrinkboxFromBand(ptr, pabyRedBand, pabyGreenBand, pabyBlueBand,
1104 : nPixels);
1105 : }
1106 5721 : else if (psHashHistogram != nullptr)
1107 : {
1108 1667 : shrinkboxFromHashHistogram(ptr, psHashHistogram);
1109 : }
1110 : else
1111 : {
1112 4054 : shrinkbox(ptr, histogram, nCLevels);
1113 : }
1114 6139 : }
1115 :
1116 : /************************************************************************/
1117 : /* shrinkbox() */
1118 : /************************************************************************/
1119 : template <class T>
1120 8120 : static void shrinkbox(Colorbox *box, const T *histogram, int nCLevels)
1121 : {
1122 8120 : if (box->rmax > box->rmin)
1123 : {
1124 5635 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1125 : {
1126 17401 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1127 : {
1128 : const T *histp =
1129 16604 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1130 222489 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1131 : {
1132 210723 : if (*histp++ != 0)
1133 : {
1134 4838 : box->rmin = ir;
1135 4838 : goto have_rmin;
1136 : }
1137 : }
1138 : }
1139 : }
1140 : }
1141 3282 : have_rmin:
1142 8120 : if (box->rmax > box->rmin)
1143 : {
1144 6338 : for (int ir = box->rmax; ir >= box->rmin; --ir)
1145 : {
1146 28343 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1147 : {
1148 : const T *histp =
1149 26815 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1150 401263 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1151 : {
1152 379258 : if (*histp++ != 0)
1153 : {
1154 4810 : box->rmax = ir;
1155 4810 : goto have_rmax;
1156 : }
1157 : }
1158 : }
1159 : }
1160 : }
1161 :
1162 3310 : have_rmax:
1163 8120 : if (box->gmax > box->gmin)
1164 : {
1165 6954 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1166 : {
1167 23477 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1168 : {
1169 : const T *histp =
1170 21999 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1171 335681 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1172 : {
1173 319156 : if (*histp++ != 0)
1174 : {
1175 5474 : box->gmin = ig;
1176 5474 : goto have_gmin;
1177 : }
1178 : }
1179 : }
1180 : }
1181 : }
1182 :
1183 2644 : have_gmin:
1184 8120 : if (box->gmax > box->gmin)
1185 : {
1186 8194 : for (int ig = box->gmax; ig >= box->gmin; --ig)
1187 : {
1188 40644 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1189 : {
1190 : const T *histp =
1191 37879 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1192 551011 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1193 : {
1194 518559 : if (*histp++ != 0)
1195 : {
1196 5427 : box->gmax = ig;
1197 5427 : goto have_gmax;
1198 : }
1199 : }
1200 : }
1201 : }
1202 : }
1203 :
1204 2691 : have_gmax:
1205 8120 : if (box->bmax > box->bmin)
1206 : {
1207 5658 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1208 : {
1209 13243 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1210 : {
1211 : const T *histp =
1212 12547 : HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
1213 129829 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1214 : {
1215 122242 : if (*histp != 0)
1216 : {
1217 4960 : box->bmin = ib;
1218 4960 : goto have_bmin;
1219 : }
1220 117282 : histp += nCLevels;
1221 : }
1222 : }
1223 : }
1224 : }
1225 :
1226 3158 : have_bmin:
1227 8120 : if (box->bmax > box->bmin)
1228 : {
1229 7414 : for (int ib = box->bmax; ib >= box->bmin; --ib)
1230 : {
1231 30629 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1232 : {
1233 : const T *histp =
1234 28139 : HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
1235 305050 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1236 : {
1237 281833 : if (*histp != 0)
1238 : {
1239 4922 : box->bmax = ib;
1240 4922 : goto have_bmax;
1241 : }
1242 276911 : histp += nCLevels;
1243 : }
1244 : }
1245 : }
1246 : }
1247 :
1248 3196 : have_bmax:;
1249 8120 : }
1250 :
1251 : // Explicitly instantiate template functions.
1252 : template int GDALComputeMedianCutPCTInternal<GUInt32>(
1253 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1254 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1255 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1256 : GUInt32 *panHistogram, GDALColorTableH hColorTable,
1257 : GDALProgressFunc pfnProgress, void *pProgressArg);
1258 :
1259 : template int GDALComputeMedianCutPCTInternal<GUIntBig>(
1260 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1261 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1262 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1263 : GUIntBig *panHistogram, GDALColorTableH hColorTable,
1264 : GDALProgressFunc pfnProgress, void *pProgressArg);
|
