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