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 18070300 : static int MAKE_COLOR_CODE(int r, int g, int b)
51 : {
52 18070300 : 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 520716 : static inline int *FindAndInsertColorCount(HashHistogram *psHashHistogram,
218 : GUInt32 nColorCode)
219 : {
220 520716 : GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
221 : while (true)
222 : {
223 520716 : if (psHashHistogram[nIdx].nColorCode == nColorCode)
224 : {
225 461027 : return &(psHashHistogram[nIdx].nCount);
226 : }
227 59689 : if (!IsColorCodeSet(psHashHistogram[nIdx].nColorCode))
228 : {
229 56308 : psHashHistogram[nIdx].nColorCode = nColorCode;
230 56308 : psHashHistogram[nIdx].nCount = 0;
231 56308 : 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 54 : 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 54 : VALIDATE_POINTER1(hRed, "GDALComputeMedianCutPCT", CE_Failure);
279 54 : VALIDATE_POINTER1(hGreen, "GDALComputeMedianCutPCT", CE_Failure);
280 54 : VALIDATE_POINTER1(hBlue, "GDALComputeMedianCutPCT", CE_Failure);
281 :
282 54 : CPLErr err = CE_None;
283 :
284 : /* -------------------------------------------------------------------- */
285 : /* Validate parameters. */
286 : /* -------------------------------------------------------------------- */
287 54 : const int nXSize = GDALGetRasterBandXSize(hRed);
288 54 : const int nYSize = GDALGetRasterBandYSize(hRed);
289 :
290 54 : if (GDALGetRasterBandXSize(hGreen) != nXSize ||
291 54 : GDALGetRasterBandYSize(hGreen) != nYSize ||
292 162 : GDALGetRasterBandXSize(hBlue) != nXSize ||
293 54 : 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 54 : 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 54 : 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 54 : 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 54 : if (pfnProgress == nullptr)
329 27 : pfnProgress = GDALDummyProgress;
330 :
331 54 : int nSrcNoData = -1;
332 54 : GDALRasterBandH hMaskBand = nullptr;
333 54 : int bSrcHasNoDataR = FALSE;
334 54 : const double dfSrcNoDataR = GDALGetRasterNoDataValue(hRed, &bSrcHasNoDataR);
335 54 : int bSrcHasNoDataG = FALSE;
336 : const double dfSrcNoDataG =
337 54 : GDALGetRasterNoDataValue(hGreen, &bSrcHasNoDataG);
338 54 : int bSrcHasNoDataB = FALSE;
339 : const double dfSrcNoDataB =
340 54 : GDALGetRasterNoDataValue(hBlue, &bSrcHasNoDataB);
341 54 : 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 52 : const int nMaskFlags = GDALGetMaskFlags(hRed);
351 52 : if ((nMaskFlags & GMF_PER_DATASET))
352 5 : hMaskBand = GDALGetMaskBand(hRed);
353 : }
354 :
355 : /* ==================================================================== */
356 : /* STEP 1: create empty boxes. */
357 : /* ==================================================================== */
358 75 : if (static_cast<GUInt32>(nXSize) >
359 54 : 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 54 : T nPixels = 0;
367 13 : if (nBits == 8 && pabyRedBand != nullptr && pabyGreenBand != nullptr &&
368 67 : pabyBlueBand != nullptr &&
369 12 : static_cast<GUInt32>(nXSize) <=
370 12 : cpl::NumericLimits<T>::max() / static_cast<GUInt32>(nYSize))
371 : {
372 12 : nPixels = static_cast<T>(nXSize) * static_cast<T>(nYSize);
373 : }
374 :
375 54 : const int nCLevels = 1 << nBits;
376 54 : const int nCLevelsCube = nCLevels * nCLevels * nCLevels;
377 54 : T *histogram = nullptr;
378 54 : HashHistogram *psHashHistogram = nullptr;
379 54 : if (panHistogram)
380 : {
381 12 : 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 11 : histogram = nullptr;
387 11 : psHashHistogram = reinterpret_cast<HashHistogram *>(panHistogram);
388 11 : 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 42 : static_cast<T *>(VSI_CALLOC_VERBOSE(nCLevelsCube, sizeof(T)));
401 42 : if (histogram == nullptr)
402 : {
403 0 : return CE_Failure;
404 : }
405 : }
406 : Colorbox *box_list =
407 54 : static_cast<Colorbox *>(CPLMalloc(nColors * sizeof(Colorbox)));
408 54 : Colorbox *freeboxes = box_list;
409 54 : freeboxes[0].next = &freeboxes[1];
410 54 : freeboxes[0].prev = nullptr;
411 12357 : for (int i = 1; i < nColors - 1; ++i)
412 : {
413 12303 : freeboxes[i].next = &freeboxes[i + 1];
414 12303 : freeboxes[i].prev = &freeboxes[i - 1];
415 : }
416 54 : freeboxes[nColors - 1].next = nullptr;
417 54 : freeboxes[nColors - 1].prev = &freeboxes[nColors - 2];
418 :
419 : /* ==================================================================== */
420 : /* Build histogram. */
421 : /* ==================================================================== */
422 :
423 : /* -------------------------------------------------------------------- */
424 : /* Initialize the box datastructures. */
425 : /* -------------------------------------------------------------------- */
426 54 : Colorbox *freeboxes_before = freeboxes;
427 54 : freeboxes = freeboxes_before->next;
428 54 : if (freeboxes)
429 54 : freeboxes->prev = nullptr;
430 :
431 54 : Colorbox *usedboxes = freeboxes_before;
432 54 : usedboxes->next = nullptr;
433 54 : usedboxes->rmin = 999;
434 54 : usedboxes->gmin = 999;
435 54 : usedboxes->bmin = 999;
436 54 : usedboxes->rmax = -1;
437 54 : usedboxes->gmax = -1;
438 54 : usedboxes->bmax = -1;
439 54 : usedboxes->total =
440 54 : 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 54 : const int nColorShift = 8 - nBits;
448 54 : int nColorCounter = 0;
449 54 : GByte anRed[256] = {};
450 54 : GByte anGreen[256] = {};
451 54 : GByte anBlue[256] = {};
452 :
453 54 : GByte *pabyRedLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
454 54 : GByte *pabyGreenLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
455 54 : GByte *pabyBlueLine = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
456 0 : std::unique_ptr<GByte, VSIFreeReleaser> pabyMask;
457 54 : if (hMaskBand)
458 5 : pabyMask.reset(static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)));
459 :
460 54 : if (pabyRedLine == nullptr || pabyGreenLine == nullptr ||
461 108 : pabyBlueLine == nullptr || (hMaskBand && !pabyMask))
462 : {
463 0 : err = CE_Failure;
464 0 : goto end_and_cleanup;
465 : }
466 :
467 9067 : for (int iLine = 0; iLine < nYSize; iLine++)
468 : {
469 9013 : 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 9013 : err = GDALRasterIO(hRed, GF_Read, 0, iLine, nXSize, 1, pabyRedLine,
478 : nXSize, 1, GDT_UInt8, 0, 0);
479 9013 : if (err == CE_None)
480 9013 : err = GDALRasterIO(hGreen, GF_Read, 0, iLine, nXSize, 1,
481 : pabyGreenLine, nXSize, 1, GDT_UInt8, 0, 0);
482 9013 : if (err == CE_None)
483 9013 : err = GDALRasterIO(hBlue, GF_Read, 0, iLine, nXSize, 1,
484 : pabyBlueLine, nXSize, 1, GDT_UInt8, 0, 0);
485 9013 : if (err == CE_None && hMaskBand)
486 1689 : err = GDALRasterIO(hMaskBand, GF_Read, 0, iLine, nXSize, 1,
487 1689 : pabyMask.get(), nXSize, 1, GDT_UInt8, 0, 0);
488 9013 : if (err != CE_None)
489 0 : goto end_and_cleanup;
490 :
491 6455690 : for (int iPixel = 0; iPixel < nXSize; iPixel++)
492 : {
493 15300680 : if ((pabyRedLine[iPixel] == nSrcNoData &&
494 33554 : pabyGreenLine[iPixel] == nSrcNoData &&
495 15276940 : pabyBlueLine[iPixel] == nSrcNoData) ||
496 8796960 : (pabyMask && (pabyMask.get())[iPixel] == 0))
497 : {
498 2373739 : continue;
499 : }
500 :
501 4072940 : const int nRed = pabyRedLine[iPixel] >> nColorShift;
502 4072940 : const int nGreen = pabyGreenLine[iPixel] >> nColorShift;
503 4072940 : const int nBlue = pabyBlueLine[iPixel] >> nColorShift;
504 :
505 4072940 : usedboxes->rmin = std::min(usedboxes->rmin, nRed);
506 4072940 : usedboxes->gmin = std::min(usedboxes->gmin, nGreen);
507 4072940 : usedboxes->bmin = std::min(usedboxes->bmin, nBlue);
508 4072940 : usedboxes->rmax = std::max(usedboxes->rmax, nRed);
509 4072940 : usedboxes->gmax = std::max(usedboxes->gmax, nGreen);
510 4072940 : usedboxes->bmax = std::max(usedboxes->bmax, nBlue);
511 :
512 : bool bFirstOccurrence;
513 4072940 : if (psHashHistogram)
514 : {
515 520716 : int *pnColor = FindAndInsertColorCount(
516 520716 : psHashHistogram, MAKE_COLOR_CODE(nRed, nGreen, nBlue));
517 520716 : bFirstOccurrence = (*pnColor == 0);
518 520716 : (*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 4072940 : if (bFirstOccurrence)
528 : {
529 111369 : if (nColorShift == 0 && nColorCounter < nColors)
530 : {
531 2316 : anRed[nColorCounter] = static_cast<GByte>(nRed);
532 2316 : anGreen[nColorCounter] = static_cast<GByte>(nGreen);
533 2316 : anBlue[nColorCounter] = static_cast<GByte>(nBlue);
534 : }
535 111369 : nColorCounter++;
536 : }
537 : }
538 : }
539 :
540 54 : 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 54 : if (nColorShift == 0 && nColorCounter <= nColors)
548 : {
549 : #if DEBUG_VERBOSE
550 : CPLDebug("MEDIAN_CUT", "%d colors found <= %d", nColorCounter, nColors);
551 : #endif
552 4 : if (panPixelCountPerColorTableEntry)
553 : {
554 1 : panPixelCountPerColorTableEntry->clear();
555 1 : panPixelCountPerColorTableEntry->reserve(nColors);
556 : }
557 16 : for (int iColor = 0; iColor < nColorCounter; iColor++)
558 : {
559 12 : const GDALColorEntry sEntry = {static_cast<GByte>(anRed[iColor]),
560 12 : static_cast<GByte>(anGreen[iColor]),
561 12 : static_cast<GByte>(anBlue[iColor]),
562 : 255};
563 12 : GDALSetColorEntry(hColorTable, iColor, &sEntry);
564 12 : 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 4 : 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 7417 : while (freeboxes != nullptr)
590 : {
591 7367 : auto ptr = largest_box(usedboxes);
592 7367 : if (ptr != nullptr)
593 7347 : splitbox(ptr, histogram, psHashHistogram, nCLevels, &freeboxes,
594 : &usedboxes, pabyRedBand, pabyGreenBand, pabyBlueBand,
595 : nPixels);
596 : else
597 20 : freeboxes = nullptr;
598 : }
599 :
600 : /* ==================================================================== */
601 : /* STEP 4: assign colors to all boxes */
602 : /* ==================================================================== */
603 50 : if (nColorCounter > 0)
604 : {
605 48 : Colorbox *ptr = usedboxes;
606 48 : if (panPixelCountPerColorTableEntry)
607 : {
608 19 : panPixelCountPerColorTableEntry->clear();
609 19 : panPixelCountPerColorTableEntry->reserve(nColors);
610 : }
611 7443 : for (int i = 0; ptr != nullptr; ++i, ptr = ptr->next)
612 : {
613 7395 : const GDALColorEntry sEntry = {
614 7395 : static_cast<GByte>(((ptr->rmin + ptr->rmax) << nColorShift) /
615 : 2),
616 7395 : static_cast<GByte>(((ptr->gmin + ptr->gmax) << nColorShift) /
617 : 2),
618 7395 : static_cast<GByte>(((ptr->bmin + ptr->bmax) << nColorShift) /
619 : 2),
620 : 255};
621 7395 : GDALSetColorEntry(hColorTable, i, &sEntry);
622 7395 : if (panPixelCountPerColorTableEntry)
623 797 : panPixelCountPerColorTableEntry->push_back(
624 797 : static_cast<T>(ptr->total));
625 : }
626 : }
627 :
628 2 : end_and_cleanup:
629 54 : CPLFree(pabyRedLine);
630 54 : CPLFree(pabyGreenLine);
631 54 : CPLFree(pabyBlueLine);
632 :
633 : // We're done with the boxes now.
634 54 : CPLFree(box_list);
635 54 : freeboxes = nullptr;
636 54 : usedboxes = nullptr;
637 :
638 54 : if (panHistogram == nullptr)
639 42 : CPLFree(histogram);
640 :
641 54 : return err;
642 : }
643 :
644 : /************************************************************************/
645 : /* largest_box() */
646 : /************************************************************************/
647 :
648 7367 : static Colorbox *largest_box(Colorbox *usedboxes)
649 : {
650 7367 : Colorbox *b = nullptr;
651 :
652 900340 : for (Colorbox *p = usedboxes; p != nullptr; p = p->next)
653 : {
654 892973 : if ((p->rmax > p->rmin || p->gmax > p->gmin || p->bmax > p->bmin) &&
655 732940 : (b == nullptr || p->total > b->total))
656 : {
657 66668 : b = p;
658 : }
659 : }
660 7367 : return b;
661 : }
662 :
663 743 : static void shrinkboxFromBand(Colorbox *ptr, const GByte *pabyRedBand,
664 : const GByte *pabyGreenBand,
665 : const GByte *pabyBlueBand, GUIntBig nPixels)
666 : {
667 743 : int rmin_new = 255;
668 743 : int rmax_new = 0;
669 743 : int gmin_new = 255;
670 743 : int gmax_new = 0;
671 743 : int bmin_new = 255;
672 743 : int bmax_new = 0;
673 35470300 : for (GUIntBig i = 0; i < nPixels; i++)
674 : {
675 35469500 : const int iR = pabyRedBand[i];
676 35469500 : const int iG = pabyGreenBand[i];
677 35469500 : const int iB = pabyBlueBand[i];
678 35469500 : if (iR >= ptr->rmin && iR <= ptr->rmax && iG >= ptr->gmin &&
679 5781920 : iG <= ptr->gmax && iB >= ptr->bmin && iB <= ptr->bmax)
680 : {
681 2154140 : if (iR < rmin_new)
682 3143 : rmin_new = iR;
683 2154140 : if (iR > rmax_new)
684 4274 : rmax_new = iR;
685 2154140 : if (iG < gmin_new)
686 3392 : gmin_new = iG;
687 2154140 : if (iG > gmax_new)
688 4635 : gmax_new = iG;
689 2154140 : if (iB < bmin_new)
690 4455 : bmin_new = iB;
691 2154140 : if (iB > bmax_new)
692 3314 : bmax_new = iB;
693 : }
694 : }
695 :
696 743 : CPLAssert(rmin_new >= ptr->rmin && rmin_new <= rmax_new &&
697 : rmax_new <= ptr->rmax);
698 743 : CPLAssert(gmin_new >= ptr->gmin && gmin_new <= gmax_new &&
699 : gmax_new <= ptr->gmax);
700 743 : CPLAssert(bmin_new >= ptr->bmin && bmin_new <= bmax_new &&
701 : bmax_new <= ptr->bmax);
702 743 : ptr->rmin = rmin_new;
703 743 : ptr->rmax = rmax_new;
704 743 : ptr->gmin = gmin_new;
705 743 : ptr->gmax = gmax_new;
706 743 : ptr->bmin = bmin_new;
707 743 : ptr->bmax = bmax_new;
708 743 : }
709 :
710 3415 : static void shrinkboxFromHashHistogram(Colorbox *box,
711 : const HashHistogram *psHashHistogram)
712 : {
713 3415 : if (box->rmax > box->rmin)
714 : {
715 4917 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
716 : {
717 54728 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
718 : {
719 770529 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
720 : {
721 720718 : if (FindColorCount(psHashHistogram,
722 1441440 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
723 : {
724 3313 : box->rmin = ir;
725 3313 : goto have_rmin;
726 : }
727 : }
728 : }
729 : }
730 : }
731 102 : have_rmin:
732 3415 : if (box->rmax > box->rmin)
733 : {
734 6050 : for (int ir = box->rmax; ir >= box->rmin; --ir)
735 : {
736 65302 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
737 : {
738 1301360 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
739 : {
740 1242110 : if (FindColorCount(psHashHistogram,
741 2484210 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
742 : {
743 3310 : box->rmax = ir;
744 3310 : goto have_rmax;
745 : }
746 : }
747 : }
748 : }
749 : }
750 :
751 105 : have_rmax:
752 3415 : if (box->gmax > box->gmin)
753 : {
754 7027 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
755 : {
756 71453 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
757 : {
758 1233930 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
759 : {
760 1169500 : if (FindColorCount(psHashHistogram,
761 2339010 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
762 : {
763 3252 : box->gmin = ig;
764 3252 : goto have_gmin;
765 : }
766 : }
767 : }
768 : }
769 : }
770 :
771 163 : have_gmin:
772 3415 : if (box->gmax > box->gmin)
773 : {
774 8709 : for (int ig = box->gmax; ig >= box->gmin; --ig)
775 : {
776 99428 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
777 : {
778 93961 : int ib = box->bmin;
779 1785560 : for (; ib <= box->bmax; ++ib)
780 : {
781 1694840 : if (FindColorCount(psHashHistogram,
782 3389690 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
783 : {
784 3242 : box->gmax = ig;
785 3242 : goto have_gmax;
786 : }
787 : }
788 : }
789 : }
790 : }
791 :
792 173 : have_gmax:
793 3415 : if (box->bmax > box->bmin)
794 : {
795 5004 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
796 : {
797 36216 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
798 : {
799 572282 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
800 : {
801 541070 : if (FindColorCount(psHashHistogram,
802 1082140 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
803 : {
804 3312 : box->bmin = ib;
805 3312 : goto have_bmin;
806 : }
807 : }
808 : }
809 : }
810 : }
811 :
812 103 : have_bmin:
813 3415 : if (box->bmax > box->bmin)
814 : {
815 6196 : for (int ib = box->bmax; ib >= box->bmin; --ib)
816 : {
817 48883 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
818 : {
819 675385 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
820 : {
821 632698 : if (FindColorCount(psHashHistogram,
822 1265400 : MAKE_COLOR_CODE(ir, ig, ib)) != 0)
823 : {
824 3298 : box->bmax = ib;
825 3298 : goto have_bmax;
826 : }
827 : }
828 : }
829 : }
830 : }
831 :
832 117 : have_bmax:;
833 3415 : }
834 :
835 : /************************************************************************/
836 : /* splitbox() */
837 : /************************************************************************/
838 : template <class T>
839 7347 : static void splitbox(Colorbox *ptr, const T *histogram,
840 : const HashHistogram *psHashHistogram, int nCLevels,
841 : Colorbox **pfreeboxes, Colorbox **pusedboxes,
842 : GByte *pabyRedBand, GByte *pabyGreenBand,
843 : GByte *pabyBlueBand, T nPixels)
844 : {
845 7347 : T hist2[256] = {};
846 7347 : int first = 0;
847 7347 : int last = 0;
848 :
849 : enum
850 : {
851 : RED,
852 : GREEN,
853 : BLUE
854 : } axis;
855 :
856 : // See which axis is the largest, do a histogram along that axis. Split at
857 : // median point. Contract both new boxes to fit points and return.
858 : {
859 7347 : int i = ptr->rmax - ptr->rmin;
860 7347 : if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
861 2466 : axis = RED;
862 4881 : else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
863 2385 : axis = GREEN;
864 : else
865 2496 : axis = BLUE;
866 : }
867 : // Get histogram along longest axis.
868 7347 : const GUInt32 nIters = (ptr->rmax - ptr->rmin + 1) *
869 7347 : (ptr->gmax - ptr->gmin + 1) *
870 7347 : (ptr->bmax - ptr->bmin + 1);
871 :
872 7347 : switch (axis)
873 : {
874 2466 : case RED:
875 : {
876 2466 : if (nPixels != 0 && nIters > nPixels)
877 : {
878 229 : const int rmin = ptr->rmin;
879 229 : const int rmax = ptr->rmax;
880 229 : const int gmin = ptr->gmin;
881 229 : const int gmax = ptr->gmax;
882 229 : const int bmin = ptr->bmin;
883 229 : const int bmax = ptr->bmax;
884 11014200 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
885 : {
886 11014000 : int iR = pabyRedBand[iPixel];
887 11014000 : int iG = pabyGreenBand[iPixel];
888 11014000 : int iB = pabyBlueBand[iPixel];
889 11014000 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
890 1379080 : iB >= bmin && iB <= bmax)
891 : {
892 893690 : hist2[iR]++;
893 : }
894 229 : }
895 : }
896 2237 : else if (psHashHistogram)
897 : {
898 472 : T *histp = &hist2[ptr->rmin];
899 11243 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
900 : {
901 10771 : *histp = 0;
902 226089 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
903 : {
904 3704290 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
905 : {
906 3488980 : *histp += FindColorCount(
907 3488980 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
908 : }
909 : }
910 10771 : histp++;
911 : }
912 : }
913 : else
914 : {
915 1765 : T *histp = &hist2[ptr->rmin];
916 16617 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
917 : {
918 14852 : *histp = 0;
919 188612 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
920 : {
921 : const T *iptr =
922 173760 : HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
923 3151039 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
924 2977279 : *histp += *iptr++;
925 : }
926 14852 : histp++;
927 : }
928 : }
929 2466 : first = ptr->rmin;
930 2466 : last = ptr->rmax;
931 2466 : break;
932 : }
933 2385 : case GREEN:
934 : {
935 2385 : if (nPixels != 0 && nIters > nPixels)
936 : {
937 153 : const int rmin = ptr->rmin;
938 153 : const int rmax = ptr->rmax;
939 153 : const int gmin = ptr->gmin;
940 153 : const int gmax = ptr->gmax;
941 153 : const int bmin = ptr->bmin;
942 153 : const int bmax = ptr->bmax;
943 7442490 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
944 : {
945 7442340 : const int iR = pabyRedBand[iPixel];
946 7442340 : const int iG = pabyGreenBand[iPixel];
947 7442340 : const int iB = pabyBlueBand[iPixel];
948 7442340 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
949 1575360 : iB >= bmin && iB <= bmax)
950 : {
951 948439 : hist2[iG]++;
952 : }
953 153 : }
954 : }
955 2232 : else if (psHashHistogram)
956 : {
957 476 : T *histp = &hist2[ptr->gmin];
958 12600 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
959 : {
960 12124 : *histp = 0;
961 183557 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
962 : {
963 3551710 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
964 : {
965 3380280 : *histp += FindColorCount(
966 3380280 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
967 : }
968 : }
969 12124 : histp++;
970 : }
971 : }
972 : else
973 : {
974 1756 : T *histp = &hist2[ptr->gmin];
975 13924 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
976 : {
977 12168 : *histp = 0;
978 110326 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
979 : {
980 : const T *iptr =
981 98158 : HISTOGRAM(histogram, nCLevels, ir, ig, ptr->bmin);
982 1913429 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
983 1815270 : *histp += *iptr++;
984 : }
985 12168 : histp++;
986 : }
987 : }
988 2385 : first = ptr->gmin;
989 2385 : last = ptr->gmax;
990 2385 : break;
991 : }
992 2496 : case BLUE:
993 : {
994 2496 : if (nPixels != 0 && nIters > nPixels)
995 : {
996 199 : const int rmin = ptr->rmin;
997 199 : const int rmax = ptr->rmax;
998 199 : const int gmin = ptr->gmin;
999 199 : const int gmax = ptr->gmax;
1000 199 : const int bmin = ptr->bmin;
1001 199 : const int bmax = ptr->bmax;
1002 9614150 : for (T iPixel = 0; iPixel < nPixels; iPixel++)
1003 : {
1004 9613950 : const int iR = pabyRedBand[iPixel];
1005 9613950 : const int iG = pabyGreenBand[iPixel];
1006 9613950 : const int iB = pabyBlueBand[iPixel];
1007 9613950 : if (iR >= rmin && iR <= rmax && iG >= gmin && iG <= gmax &&
1008 832271 : iB >= bmin && iB <= bmax)
1009 : {
1010 649243 : hist2[iB]++;
1011 : }
1012 199 : }
1013 : }
1014 2297 : else if (psHashHistogram)
1015 : {
1016 573 : T *histp = &hist2[ptr->bmin];
1017 16540 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
1018 : {
1019 15967 : *histp = 0;
1020 252363 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
1021 : {
1022 4915820 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
1023 : {
1024 4679430 : *histp += FindColorCount(
1025 4679430 : psHashHistogram, MAKE_COLOR_CODE(ir, ig, ib));
1026 : }
1027 : }
1028 15967 : histp++;
1029 : }
1030 : }
1031 : else
1032 : {
1033 1724 : T *histp = &hist2[ptr->bmin];
1034 15464 : for (int ib = ptr->bmin; ib <= ptr->bmax; ++ib)
1035 : {
1036 13740 : *histp = 0;
1037 136611 : for (int ir = ptr->rmin; ir <= ptr->rmax; ++ir)
1038 : {
1039 : const T *iptr =
1040 122871 : HISTOGRAM(histogram, nCLevels, ir, ptr->gmin, ib);
1041 2358380 : for (int ig = ptr->gmin; ig <= ptr->gmax; ++ig)
1042 : {
1043 2235514 : *histp += *iptr;
1044 2235514 : iptr += nCLevels;
1045 : }
1046 : }
1047 13740 : histp++;
1048 : }
1049 : }
1050 2496 : first = ptr->bmin;
1051 2496 : last = ptr->bmax;
1052 2496 : break;
1053 : }
1054 : }
1055 : // Find median point.
1056 7347 : T *histp = &hist2[first];
1057 7347 : int i = first; // TODO(schwehr): Rename i.
1058 : {
1059 7347 : T sum = 0;
1060 7347 : T sum2 = static_cast<T>(ptr->total / 2);
1061 57172 : for (; i <= last && (sum += *histp++) < sum2; ++i)
1062 : {
1063 : }
1064 : }
1065 7347 : if (i == first)
1066 1799 : i++;
1067 :
1068 : // Create new box, re-allocate points.
1069 7347 : Colorbox *new_cb = *pfreeboxes;
1070 7347 : *pfreeboxes = new_cb->next;
1071 7347 : if (*pfreeboxes)
1072 7317 : (*pfreeboxes)->prev = nullptr;
1073 7347 : if (*pusedboxes)
1074 7347 : (*pusedboxes)->prev = new_cb;
1075 7347 : new_cb->next = *pusedboxes;
1076 7347 : *pusedboxes = new_cb;
1077 :
1078 7347 : histp = &hist2[first];
1079 : {
1080 7347 : T sum1 = 0;
1081 58971 : for (int j = first; j < i; j++)
1082 51624 : sum1 += *histp++;
1083 7347 : T sum2 = 0;
1084 91652 : for (int j = i; j <= last; j++)
1085 84305 : sum2 += *histp++;
1086 7347 : new_cb->total = sum1;
1087 7347 : ptr->total = sum2;
1088 : }
1089 7347 : new_cb->rmin = ptr->rmin;
1090 7347 : new_cb->rmax = ptr->rmax;
1091 7347 : new_cb->gmin = ptr->gmin;
1092 7347 : new_cb->gmax = ptr->gmax;
1093 7347 : new_cb->bmin = ptr->bmin;
1094 7347 : new_cb->bmax = ptr->bmax;
1095 7347 : switch (axis)
1096 : {
1097 2466 : case RED:
1098 2466 : new_cb->rmax = i - 1;
1099 2466 : ptr->rmin = i;
1100 2466 : break;
1101 2385 : case GREEN:
1102 2385 : new_cb->gmax = i - 1;
1103 2385 : ptr->gmin = i;
1104 2385 : break;
1105 2496 : case BLUE:
1106 2496 : new_cb->bmax = i - 1;
1107 2496 : ptr->bmin = i;
1108 2496 : break;
1109 : }
1110 :
1111 7347 : if (nPixels != 0 &&
1112 2295 : static_cast<T>(new_cb->rmax - new_cb->rmin + 1) *
1113 2295 : static_cast<T>(new_cb->gmax - new_cb->gmin + 1) *
1114 2295 : static_cast<T>(new_cb->bmax - new_cb->bmin + 1) >
1115 : nPixels)
1116 : {
1117 325 : shrinkboxFromBand(new_cb, pabyRedBand, pabyGreenBand, pabyBlueBand,
1118 : nPixels);
1119 : }
1120 7022 : else if (psHashHistogram != nullptr)
1121 : {
1122 1748 : shrinkboxFromHashHistogram(new_cb, psHashHistogram);
1123 : }
1124 : else
1125 : {
1126 5274 : shrinkbox(new_cb, histogram, nCLevels);
1127 : }
1128 :
1129 7347 : if (nPixels != 0 && static_cast<T>(ptr->rmax - ptr->rmin + 1) *
1130 2295 : static_cast<T>(ptr->gmax - ptr->gmin + 1) *
1131 2295 : static_cast<T>(ptr->bmax - ptr->bmin + 1) >
1132 : nPixels)
1133 : {
1134 418 : shrinkboxFromBand(ptr, pabyRedBand, pabyGreenBand, pabyBlueBand,
1135 : nPixels);
1136 : }
1137 6929 : else if (psHashHistogram != nullptr)
1138 : {
1139 1667 : shrinkboxFromHashHistogram(ptr, psHashHistogram);
1140 : }
1141 : else
1142 : {
1143 5262 : shrinkbox(ptr, histogram, nCLevels);
1144 : }
1145 7347 : }
1146 :
1147 : /************************************************************************/
1148 : /* shrinkbox() */
1149 : /************************************************************************/
1150 : template <class T>
1151 10536 : static void shrinkbox(Colorbox *box, const T *histogram, int nCLevels)
1152 : {
1153 10536 : if (box->rmax > box->rmin)
1154 : {
1155 7042 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1156 : {
1157 20590 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1158 : {
1159 : const T *histp =
1160 19589 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1161 244021 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1162 : {
1163 230473 : if (*histp++ != 0)
1164 : {
1165 6041 : box->rmin = ir;
1166 6041 : goto have_rmin;
1167 : }
1168 : }
1169 : }
1170 : }
1171 : }
1172 4495 : have_rmin:
1173 10536 : if (box->rmax > box->rmin)
1174 : {
1175 7820 : for (int ir = box->rmax; ir >= box->rmin; --ir)
1176 : {
1177 34920 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1178 : {
1179 : const T *histp =
1180 33087 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1181 463269 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1182 : {
1183 436169 : if (*histp++ != 0)
1184 : {
1185 5987 : box->rmax = ir;
1186 5987 : goto have_rmax;
1187 : }
1188 : }
1189 : }
1190 : }
1191 : }
1192 :
1193 4549 : have_rmax:
1194 10536 : if (box->gmax > box->gmin)
1195 : {
1196 9144 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1197 : {
1198 29845 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1199 : {
1200 : const T *histp =
1201 27743 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1202 401339 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1203 : {
1204 380620 : if (*histp++ != 0)
1205 : {
1206 7024 : box->gmin = ig;
1207 7024 : goto have_gmin;
1208 : }
1209 : }
1210 : }
1211 : }
1212 : }
1213 :
1214 3494 : have_gmin:
1215 10536 : if (box->gmax > box->gmin)
1216 : {
1217 11393 : for (int ig = box->gmax; ig >= box->gmin; --ig)
1218 : {
1219 51896 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1220 : {
1221 : const T *histp =
1222 47405 : HISTOGRAM(histogram, nCLevels, ir, ig, box->bmin);
1223 666264 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1224 : {
1225 625743 : if (*histp++ != 0)
1226 : {
1227 6884 : box->gmax = ig;
1228 6884 : goto have_gmax;
1229 : }
1230 : }
1231 : }
1232 : }
1233 : }
1234 :
1235 3634 : have_gmax:
1236 10536 : if (box->bmax > box->bmin)
1237 : {
1238 7845 : for (int ib = box->bmin; ib <= box->bmax; ++ib)
1239 : {
1240 18296 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1241 : {
1242 : const T *histp =
1243 17082 : HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
1244 166256 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1245 : {
1246 155789 : if (*histp != 0)
1247 : {
1248 6615 : box->bmin = ib;
1249 6615 : goto have_bmin;
1250 : }
1251 149174 : histp += nCLevels;
1252 : }
1253 : }
1254 : }
1255 : }
1256 :
1257 3905 : have_bmin:
1258 10536 : if (box->bmax > box->bmin)
1259 : {
1260 10612 : for (int ib = box->bmax; ib >= box->bmin; --ib)
1261 : {
1262 40412 : for (int ir = box->rmin; ir <= box->rmax; ++ir)
1263 : {
1264 : const T *histp =
1265 36299 : HISTOGRAM(histogram, nCLevels, ir, box->gmin, ib);
1266 379858 : for (int ig = box->gmin; ig <= box->gmax; ++ig)
1267 : {
1268 350042 : if (*histp != 0)
1269 : {
1270 6483 : box->bmax = ib;
1271 6483 : goto have_bmax;
1272 : }
1273 343559 : histp += nCLevels;
1274 : }
1275 : }
1276 : }
1277 : }
1278 :
1279 4037 : have_bmax:;
1280 10536 : }
1281 :
1282 : // Explicitly instantiate template functions.
1283 : template int GDALComputeMedianCutPCTInternal<GUInt32>(
1284 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1285 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1286 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1287 : GUInt32 *panHistogram, GDALColorTableH hColorTable,
1288 : GDALProgressFunc pfnProgress, void *pProgressArg,
1289 : std::vector<GUInt32> *panPixelCountPerColorTableEntry);
1290 :
1291 : template int GDALComputeMedianCutPCTInternal<GUIntBig>(
1292 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1293 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1294 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1295 : GUIntBig *panHistogram, GDALColorTableH hColorTable,
1296 : GDALProgressFunc pfnProgress, void *pProgressArg,
1297 : std::vector<GUIntBig> *panPixelCountPerColorTableEntry);
1298 :
1299 0 : int GDALComputeMedianCutPCT(
1300 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1301 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1302 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1303 : GUInt32 *panHistogram, GDALColorTableH hColorTable,
1304 : GDALProgressFunc pfnProgress, void *pProgressArg,
1305 : std::vector<GUInt32> *panPixelCountPerColorTableEntry)
1306 : {
1307 0 : return GDALComputeMedianCutPCTInternal<GUInt32>(
1308 : hRed, hGreen, hBlue, pabyRedBand, pabyGreenBand, pabyBlueBand,
1309 : pfnIncludePixel, nColors, nBits, panHistogram, hColorTable, pfnProgress,
1310 0 : pProgressArg, panPixelCountPerColorTableEntry);
1311 : }
1312 :
1313 21 : int GDALComputeMedianCutPCT(
1314 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
1315 : GByte *pabyRedBand, GByte *pabyGreenBand, GByte *pabyBlueBand,
1316 : int (*pfnIncludePixel)(int, int, void *), int nColors, int nBits,
1317 : GUIntBig *panHistogram, GDALColorTableH hColorTable,
1318 : GDALProgressFunc pfnProgress, void *pProgressArg,
1319 : std::vector<GUIntBig> *panPixelCountPerColorTableEntry)
1320 : {
1321 21 : return GDALComputeMedianCutPCTInternal<GUIntBig>(
1322 : hRed, hGreen, hBlue, pabyRedBand, pabyGreenBand, pabyBlueBand,
1323 : pfnIncludePixel, nColors, nBits, panHistogram, hColorTable, pfnProgress,
1324 21 : pProgressArg, panPixelCountPerColorTableEntry);
1325 : }
|
