Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: CIETMap Phase 2
4 : * Purpose: Convert RGB (24bit) to a pseudo-colored approximation using
5 : * Floyd-Steinberg dithering (error diffusion).
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2001, Frank Warmerdam
10 : * Copyright (c) 2007, Even Rouault <even dot rouault at spatialys.com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ******************************************************************************
14 : *
15 : * Notes:
16 : *
17 : * [1] Floyd-Steinberg dither:
18 : * I should point out that the actual fractions we used were, assuming
19 : * you are at X, moving left to right:
20 : *
21 : * X 7/16
22 : * 3/16 5/16 1/16
23 : *
24 : * Note that the error goes to four neighbors, not three. I think this
25 : * will probably do better (at least for black and white) than the
26 : * 3/8-3/8-1/4 distribution, at the cost of greater processing. I have
27 : * seen the 3/8-3/8-1/4 distribution described as "our" algorithm before,
28 : * but I have no idea who the credit really belongs to.
29 : * --
30 : * Lou Steinberg
31 : */
32 :
33 : #include "cpl_port.h"
34 : #include "gdal_alg.h"
35 : #include "gdal_alg_priv.h"
36 :
37 : #include <cmath>
38 : #include <cstdlib>
39 : #include <cstring>
40 : #include <algorithm>
41 :
42 : #include "cpl_conv.h"
43 : #include "cpl_error.h"
44 : #include "cpl_progress.h"
45 : #include "cpl_vsi.h"
46 : #include "gdal.h"
47 : #include "gdal_priv.h"
48 :
49 : #ifdef USE_NEON_OPTIMIZATIONS
50 : #define USE_SSE2
51 : #include "include_sse2neon.h"
52 : #elif defined(__x86_64) || defined(_M_X64)
53 : #define USE_SSE2
54 : #include <emmintrin.h>
55 : #endif
56 :
57 : #ifdef USE_SSE2
58 :
59 : #define CAST_PCT(x) reinterpret_cast<GByte *>(x)
60 : #define ALIGN_INT_ARRAY_ON_16_BYTE(x) \
61 : (((reinterpret_cast<GUIntptr_t>(x) % 16) != 0) \
62 : ? reinterpret_cast<int *>(reinterpret_cast<GByte *>(x) + 16 - \
63 : (reinterpret_cast<GUIntptr_t>(x) % 16)) \
64 : : (x))
65 : #else
66 : #define CAST_PCT(x) x
67 : #endif
68 :
69 : #ifndef MAKE_COLOR_CODE_defined
70 : #define MAKE_COLOR_CODE_defined
71 :
72 1345150 : static int MAKE_COLOR_CODE(int r, int g, int b)
73 : {
74 1345150 : return r | (g << 8) | (b << 16);
75 : }
76 : #endif
77 :
78 : static void FindNearestColor(int nColors, const int *panPCT,
79 : GByte *pabyColorMap, int nCLevels, int nDstNoData);
80 : static int FindNearestColor(int nColors, const int *panPCT, int nRedValue,
81 : int nGreenValue, int nBlueValue, int nDstNoData);
82 :
83 : // Structure for a hashmap from a color code to a color index of the
84 : // color table.
85 :
86 : // NOTE: if changing the size of this structure, edit
87 : // MEDIAN_CUT_AND_DITHER_BUFFER_SIZE_65536 in gdal_alg_priv.h and take
88 : // into account HashHistogram in gdalmediancut.cpp.
89 : typedef struct
90 : {
91 : GUInt32 nColorCode;
92 : GUInt32 nColorCode2;
93 : GUInt32 nColorCode3;
94 : GByte nIndex;
95 : GByte nIndex2;
96 : GByte nIndex3;
97 : GByte nPadding;
98 : } ColorIndex;
99 :
100 : /************************************************************************/
101 : /* GDALDitherRGB2PCT() */
102 : /************************************************************************/
103 :
104 : #ifndef IsColorCodeSet_defined
105 : #define IsColorCodeSet_defined
106 :
107 146386 : static inline bool IsColorCodeSet(GUInt32 nColorCode)
108 : {
109 146386 : return (nColorCode >> 31) == 0;
110 : }
111 : #endif
112 :
113 : /**
114 : * 24bit to 8bit conversion with dithering.
115 : *
116 : * This functions utilizes Floyd-Steinberg dithering in the process of
117 : * converting a 24bit RGB image into a pseudocolored 8bit image using a
118 : * provided color table.
119 : *
120 : * The red, green and blue input bands do not necessarily need to come
121 : * from the same file, but they must be the same width and height. They will
122 : * be clipped to 8bit during reading, so non-eight bit bands are generally
123 : * inappropriate. Likewise the hTarget band will be written with 8bit values
124 : * and must match the width and height of the source bands.
125 : *
126 : * The color table cannot have more than 256 entries.
127 : *
128 : * Since GDAL 3.13, source nodata values or mask band will be taken into account
129 : * to determine which pixels are valid, provided that a destination nodata value
130 : * is set to hTarget.
131 : *
132 : * @param hRed Red input band.
133 : * @param hGreen Green input band.
134 : * @param hBlue Blue input band.
135 : * @param hTarget Output band.
136 : * @param hColorTable the color table to use with the output band.
137 : * @param pfnProgress callback for reporting algorithm progress matching the
138 : * GDALProgressFunc() semantics. May be NULL.
139 : * @param pProgressArg callback argument passed to pfnProgress.
140 : *
141 : * @return CE_None on success or CE_Failure if an error occurs.
142 : */
143 :
144 7 : int CPL_STDCALL GDALDitherRGB2PCT(GDALRasterBandH hRed, GDALRasterBandH hGreen,
145 : GDALRasterBandH hBlue,
146 : GDALRasterBandH hTarget,
147 : GDALColorTableH hColorTable,
148 : GDALProgressFunc pfnProgress,
149 : void *pProgressArg)
150 :
151 : {
152 7 : return GDALDitherRGB2PCTInternal(hRed, hGreen, hBlue, hTarget, hColorTable,
153 : 5, nullptr, TRUE, pfnProgress,
154 7 : pProgressArg);
155 : }
156 :
157 32 : int GDALDitherRGB2PCTInternal(
158 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
159 : GDALRasterBandH hTarget, GDALColorTableH hColorTable, int nBits,
160 : // NULL or at least 256 * 256 * 256 * sizeof(GInt16) bytes.
161 : GInt16 *pasDynamicColorMap, int bDither, GDALProgressFunc pfnProgress,
162 : void *pProgressArg)
163 : {
164 32 : VALIDATE_POINTER1(hRed, "GDALDitherRGB2PCT", CE_Failure);
165 32 : VALIDATE_POINTER1(hGreen, "GDALDitherRGB2PCT", CE_Failure);
166 32 : VALIDATE_POINTER1(hBlue, "GDALDitherRGB2PCT", CE_Failure);
167 32 : VALIDATE_POINTER1(hTarget, "GDALDitherRGB2PCT", CE_Failure);
168 32 : VALIDATE_POINTER1(hColorTable, "GDALDitherRGB2PCT", CE_Failure);
169 :
170 : /* -------------------------------------------------------------------- */
171 : /* Validate parameters. */
172 : /* -------------------------------------------------------------------- */
173 32 : const int nXSize = GDALGetRasterBandXSize(hRed);
174 32 : const int nYSize = GDALGetRasterBandYSize(hRed);
175 :
176 32 : if (GDALGetRasterBandXSize(hGreen) != nXSize ||
177 32 : GDALGetRasterBandYSize(hGreen) != nYSize ||
178 96 : GDALGetRasterBandXSize(hBlue) != nXSize ||
179 32 : GDALGetRasterBandYSize(hBlue) != nYSize)
180 : {
181 0 : CPLError(CE_Failure, CPLE_IllegalArg,
182 : "Green or blue band doesn't match size of red band.");
183 :
184 0 : return CE_Failure;
185 : }
186 :
187 64 : if (GDALGetRasterBandXSize(hTarget) != nXSize ||
188 32 : GDALGetRasterBandYSize(hTarget) != nYSize)
189 : {
190 0 : CPLError(CE_Failure, CPLE_IllegalArg,
191 : "GDALDitherRGB2PCT(): "
192 : "Target band doesn't match size of source bands.");
193 :
194 0 : return CE_Failure;
195 : }
196 :
197 32 : if (pfnProgress == nullptr)
198 23 : pfnProgress = GDALDummyProgress;
199 :
200 32 : int nSrcNoData = -1;
201 32 : GDALRasterBandH hMaskBand = nullptr;
202 32 : GByte byDstNoData = 0;
203 32 : int nDstNoData = -1;
204 32 : int bDstHasNoData = FALSE;
205 : const double dfDstNoData =
206 32 : GDALGetRasterNoDataValue(hTarget, &bDstHasNoData);
207 32 : if (bDstHasNoData && dfDstNoData >= 0 && dfDstNoData <= 255 &&
208 3 : std::round(dfDstNoData) == dfDstNoData)
209 : {
210 3 : byDstNoData = static_cast<GByte>(dfDstNoData);
211 3 : nDstNoData = byDstNoData;
212 :
213 3 : int bSrcHasNoDataR = FALSE;
214 : const double dfSrcNoDataR =
215 3 : GDALGetRasterNoDataValue(hRed, &bSrcHasNoDataR);
216 3 : int bSrcHasNoDataG = FALSE;
217 : const double dfSrcNoDataG =
218 3 : GDALGetRasterNoDataValue(hGreen, &bSrcHasNoDataG);
219 3 : int bSrcHasNoDataB = FALSE;
220 : const double dfSrcNoDataB =
221 3 : GDALGetRasterNoDataValue(hBlue, &bSrcHasNoDataB);
222 3 : if (bSrcHasNoDataR && bSrcHasNoDataG && bSrcHasNoDataB &&
223 2 : dfSrcNoDataR == dfSrcNoDataG && dfSrcNoDataR == dfSrcNoDataB &&
224 2 : dfSrcNoDataR >= 0 && dfSrcNoDataR <= 255 &&
225 2 : std::round(dfSrcNoDataR) == dfSrcNoDataR)
226 : {
227 2 : nSrcNoData = static_cast<int>(dfSrcNoDataR);
228 : }
229 : else
230 : {
231 1 : const int nMaskFlags = GDALGetMaskFlags(hRed);
232 1 : if ((nMaskFlags & GMF_PER_DATASET))
233 1 : hMaskBand = GDALGetMaskBand(hRed);
234 : }
235 : }
236 :
237 : /* -------------------------------------------------------------------- */
238 : /* Setup more direct colormap. */
239 : /* -------------------------------------------------------------------- */
240 : int iColor;
241 : #ifdef USE_SSE2
242 : int anPCTUnaligned[256 + 4]; // 4 for alignment on 16-byte boundary.
243 32 : int *anPCT = ALIGN_INT_ARRAY_ON_16_BYTE(anPCTUnaligned);
244 : #else
245 : int anPCT[256 * 4] = {};
246 : #endif
247 32 : const int nColors = GDALGetColorEntryCount(hColorTable);
248 :
249 32 : if (nColors == 0)
250 : {
251 0 : CPLError(CE_Failure, CPLE_IllegalArg,
252 : "GDALDitherRGB2PCT(): "
253 : "Color table must not be empty.");
254 :
255 0 : return CE_Failure;
256 : }
257 32 : else if (nColors > 256)
258 : {
259 0 : CPLError(CE_Failure, CPLE_IllegalArg,
260 : "GDALDitherRGB2PCT(): "
261 : "Color table cannot have more than 256 entries.");
262 :
263 0 : return CE_Failure;
264 : }
265 :
266 32 : iColor = 0;
267 7417 : do
268 : {
269 : GDALColorEntry sEntry;
270 :
271 7449 : GDALGetColorEntryAsRGB(hColorTable, iColor, &sEntry);
272 7449 : CAST_PCT(anPCT)[4 * iColor + 0] = static_cast<GByte>(sEntry.c1);
273 7449 : CAST_PCT(anPCT)[4 * iColor + 1] = static_cast<GByte>(sEntry.c2);
274 7449 : CAST_PCT(anPCT)[4 * iColor + 2] = static_cast<GByte>(sEntry.c3);
275 7449 : CAST_PCT(anPCT)[4 * iColor + 3] = 0;
276 :
277 7449 : iColor++;
278 7449 : } while (iColor < nColors);
279 :
280 : #ifdef USE_SSE2
281 : // Pad to multiple of 8 colors.
282 32 : const int nColorsMod8 = nColors % 8;
283 32 : if (nColorsMod8)
284 : {
285 1 : int iDest = nColors;
286 8 : for (iColor = 0; iColor < 8 - nColorsMod8 && iDest < 256;
287 : iColor++, iDest++)
288 : {
289 7 : anPCT[iDest] = anPCT[nColors - 1];
290 : }
291 : }
292 : #endif
293 :
294 : /* -------------------------------------------------------------------- */
295 : /* Setup various variables. */
296 : /* -------------------------------------------------------------------- */
297 32 : const int nCLevels = 1 << nBits;
298 32 : const int nCLevelsCube = nCLevels * nCLevels * nCLevels;
299 32 : ColorIndex *psColorIndexMap = nullptr;
300 :
301 32 : GByte *pabyRed = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
302 32 : GByte *pabyGreen = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
303 32 : GByte *pabyBlue = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
304 32 : std::unique_ptr<GByte, VSIFreeReleaser> pabyMask;
305 32 : if (hMaskBand)
306 1 : pabyMask.reset(static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)));
307 :
308 32 : GByte *pabyIndex = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
309 :
310 : int *panError =
311 32 : static_cast<int *>(VSI_CALLOC_VERBOSE(sizeof(int), (nXSize + 2) * 3));
312 :
313 32 : if (pabyRed == nullptr || pabyGreen == nullptr || pabyBlue == nullptr ||
314 64 : pabyIndex == nullptr || panError == nullptr || (hMaskBand && !pabyMask))
315 : {
316 0 : CPLFree(pabyRed);
317 0 : CPLFree(pabyGreen);
318 0 : CPLFree(pabyBlue);
319 0 : CPLFree(pabyIndex);
320 0 : CPLFree(panError);
321 :
322 0 : return CE_Failure;
323 : }
324 :
325 32 : GByte *pabyColorMap = nullptr;
326 32 : if (pasDynamicColorMap == nullptr)
327 : {
328 : /* --------------------------------------------------------------------
329 : */
330 : /* Build a 24bit to 8 bit color mapping. */
331 : /* --------------------------------------------------------------------
332 : */
333 :
334 : pabyColorMap = static_cast<GByte *>(
335 22 : VSI_MALLOC_VERBOSE(nCLevelsCube * sizeof(GByte)));
336 22 : if (pabyColorMap == nullptr)
337 : {
338 0 : CPLFree(pabyRed);
339 0 : CPLFree(pabyGreen);
340 0 : CPLFree(pabyBlue);
341 0 : CPLFree(pabyIndex);
342 0 : CPLFree(panError);
343 0 : CPLFree(pabyColorMap);
344 :
345 0 : return CE_Failure;
346 : }
347 :
348 22 : FindNearestColor(nColors, anPCT, pabyColorMap, nCLevels, nDstNoData);
349 : }
350 : else
351 : {
352 10 : pabyColorMap = nullptr;
353 10 : if (nBits == 8 && static_cast<GIntBig>(nXSize) * nYSize <= 65536)
354 : {
355 : // If the image is small enough, then the number of colors
356 : // will be limited and using a hashmap, rather than a full table
357 : // will be more efficient.
358 9 : psColorIndexMap =
359 : reinterpret_cast<ColorIndex *>(pasDynamicColorMap);
360 9 : memset(psColorIndexMap, 0xFF, sizeof(ColorIndex) * PRIME_FOR_65536);
361 : }
362 : else
363 : {
364 1 : memset(pasDynamicColorMap, 0xFF, 256 * 256 * 256 * sizeof(GInt16));
365 : }
366 : }
367 :
368 : /* ==================================================================== */
369 : /* Loop over all scanlines of data to process. */
370 : /* ==================================================================== */
371 32 : CPLErr err = CE_None;
372 :
373 5206 : for (int iScanline = 0; iScanline < nYSize; iScanline++)
374 : {
375 : /* --------------------------------------------------------------------
376 : */
377 : /* Report progress */
378 : /* --------------------------------------------------------------------
379 : */
380 5174 : if (!pfnProgress(iScanline / static_cast<double>(nYSize), nullptr,
381 : pProgressArg))
382 : {
383 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User Terminated");
384 0 : CPLFree(pabyRed);
385 0 : CPLFree(pabyGreen);
386 0 : CPLFree(pabyBlue);
387 0 : CPLFree(pabyIndex);
388 0 : CPLFree(panError);
389 0 : CPLFree(pabyColorMap);
390 :
391 0 : return CE_Failure;
392 : }
393 :
394 : /* --------------------------------------------------------------------
395 : */
396 : /* Read source data. */
397 : /* --------------------------------------------------------------------
398 : */
399 5174 : CPLErr err1 = GDALRasterIO(hRed, GF_Read, 0, iScanline, nXSize, 1,
400 : pabyRed, nXSize, 1, GDT_UInt8, 0, 0);
401 5174 : if (err1 == CE_None)
402 5174 : err1 = GDALRasterIO(hGreen, GF_Read, 0, iScanline, nXSize, 1,
403 : pabyGreen, nXSize, 1, GDT_UInt8, 0, 0);
404 5174 : if (err1 == CE_None)
405 5174 : err1 = GDALRasterIO(hBlue, GF_Read, 0, iScanline, nXSize, 1,
406 : pabyBlue, nXSize, 1, GDT_UInt8, 0, 0);
407 5174 : if (err1 == CE_None && hMaskBand)
408 150 : err1 = GDALRasterIO(hMaskBand, GF_Read, 0, iScanline, nXSize, 1,
409 150 : pabyMask.get(), nXSize, 1, GDT_UInt8, 0, 0);
410 5174 : if (err1 != CE_None)
411 : {
412 0 : CPLFree(pabyRed);
413 0 : CPLFree(pabyGreen);
414 0 : CPLFree(pabyBlue);
415 0 : CPLFree(pabyIndex);
416 0 : CPLFree(panError);
417 0 : CPLFree(pabyColorMap);
418 :
419 0 : return err1;
420 : }
421 :
422 : /* --------------------------------------------------------------------
423 : */
424 : /* Apply the error from the previous line to this one. */
425 : /* --------------------------------------------------------------------
426 : */
427 5174 : if (bDither)
428 : {
429 1236570 : for (int i = 0; i < nXSize; i++)
430 : {
431 1305440 : if (nDstNoData >= 0 &&
432 72900 : ((pabyRed[i] == nSrcNoData && pabyGreen[i] == nSrcNoData &&
433 72900 : pabyBlue[i] == nSrcNoData) ||
434 63892 : (pabyMask && (pabyMask.get())[i] == 0)))
435 : {
436 45088 : continue;
437 : }
438 :
439 1187460 : pabyRed[i] = static_cast<GByte>(std::max(
440 1187460 : 0, std::min(255, (pabyRed[i] + panError[i * 3 + 0 + 3]))));
441 1187460 : pabyGreen[i] = static_cast<GByte>(std::max(
442 2374910 : 0,
443 1187460 : std::min(255, (pabyGreen[i] + panError[i * 3 + 1 + 3]))));
444 1187460 : pabyBlue[i] = static_cast<GByte>(std::max(
445 1187460 : 0, std::min(255, (pabyBlue[i] + panError[i * 3 + 2 + 3]))));
446 : }
447 :
448 4024 : memset(panError, 0, sizeof(int) * (nXSize + 2) * 3);
449 : }
450 :
451 : /* --------------------------------------------------------------------
452 : */
453 : /* Figure out the nearest color to the RGB value. */
454 : /* --------------------------------------------------------------------
455 : */
456 5174 : int nLastRedError = 0;
457 5174 : int nLastGreenError = 0;
458 5174 : int nLastBlueError = 0;
459 :
460 1525220 : for (int i = 0; i < nXSize; i++)
461 : {
462 1592940 : if (nDstNoData >= 0 &&
463 72900 : ((pabyRed[i] == nSrcNoData && pabyGreen[i] == nSrcNoData &&
464 72902 : pabyBlue[i] == nSrcNoData) ||
465 63892 : (pabyMask && (pabyMask.get())[i] == 0)))
466 : {
467 45088 : nLastRedError = 0;
468 45088 : nLastGreenError = 0;
469 45088 : nLastBlueError = 0;
470 45088 : pabyIndex[i] = byDstNoData;
471 45088 : continue;
472 : }
473 :
474 : const int nRedValue =
475 1474960 : std::max(0, std::min(255, pabyRed[i] + nLastRedError));
476 : const int nGreenValue =
477 1474960 : std::max(0, std::min(255, pabyGreen[i] + nLastGreenError));
478 : const int nBlueValue =
479 1474960 : std::max(0, std::min(255, pabyBlue[i] + nLastBlueError));
480 :
481 1474960 : int iIndex = 0;
482 1474960 : int nError = 0;
483 1474960 : int nSixth = 0;
484 1474960 : if (psColorIndexMap)
485 : {
486 : const GUInt32 nColorCode =
487 389644 : MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue);
488 389644 : GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
489 : while (true)
490 : {
491 389651 : if (psColorIndexMap[nIdx].nColorCode == nColorCode)
492 : {
493 252924 : iIndex = psColorIndexMap[nIdx].nIndex;
494 252924 : break;
495 : }
496 136727 : if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode))
497 : {
498 125973 : psColorIndexMap[nIdx].nColorCode = nColorCode;
499 125973 : iIndex = FindNearestColor(nColors, anPCT, nRedValue,
500 : nGreenValue, nBlueValue,
501 : nDstNoData);
502 125973 : psColorIndexMap[nIdx].nIndex =
503 : static_cast<GByte>(iIndex);
504 125973 : break;
505 : }
506 10754 : if (psColorIndexMap[nIdx].nColorCode2 == nColorCode)
507 : {
508 1476 : iIndex = psColorIndexMap[nIdx].nIndex2;
509 1476 : break;
510 : }
511 9278 : if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode2))
512 : {
513 8876 : psColorIndexMap[nIdx].nColorCode2 = nColorCode;
514 8876 : iIndex = FindNearestColor(nColors, anPCT, nRedValue,
515 : nGreenValue, nBlueValue,
516 : nDstNoData);
517 8876 : psColorIndexMap[nIdx].nIndex2 =
518 : static_cast<GByte>(iIndex);
519 8876 : break;
520 : }
521 402 : if (psColorIndexMap[nIdx].nColorCode3 == nColorCode)
522 : {
523 31 : iIndex = psColorIndexMap[nIdx].nIndex3;
524 31 : break;
525 : }
526 371 : if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode3))
527 : {
528 364 : psColorIndexMap[nIdx].nColorCode3 = nColorCode;
529 364 : iIndex = FindNearestColor(nColors, anPCT, nRedValue,
530 : nGreenValue, nBlueValue,
531 : nDstNoData);
532 364 : psColorIndexMap[nIdx].nIndex3 =
533 : static_cast<GByte>(iIndex);
534 364 : break;
535 : }
536 :
537 0 : do
538 : {
539 7 : nIdx += 257;
540 7 : if (nIdx >= PRIME_FOR_65536)
541 0 : nIdx -= PRIME_FOR_65536;
542 : } while (
543 7 : IsColorCodeSet(psColorIndexMap[nIdx].nColorCode) &&
544 6 : psColorIndexMap[nIdx].nColorCode != nColorCode &&
545 3 : IsColorCodeSet(psColorIndexMap[nIdx].nColorCode2) &&
546 0 : psColorIndexMap[nIdx].nColorCode2 != nColorCode &&
547 10 : IsColorCodeSet(psColorIndexMap[nIdx].nColorCode3) &&
548 0 : psColorIndexMap[nIdx].nColorCode3 != nColorCode);
549 : }
550 : }
551 1085310 : else if (pasDynamicColorMap == nullptr)
552 : {
553 1005310 : const int iRed = nRedValue * nCLevels / 256;
554 1005310 : const int iGreen = nGreenValue * nCLevels / 256;
555 1005310 : const int iBlue = nBlueValue * nCLevels / 256;
556 :
557 1005310 : iIndex = pabyColorMap[iRed + iGreen * nCLevels +
558 1005310 : iBlue * nCLevels * nCLevels];
559 : }
560 : else
561 : {
562 : const GUInt32 nColorCode =
563 80000 : MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue);
564 80000 : GInt16 *psIndex = &pasDynamicColorMap[nColorCode];
565 80000 : if (*psIndex < 0)
566 : {
567 19399 : *psIndex = static_cast<GInt16>(
568 19399 : FindNearestColor(nColors, anPCT, nRedValue, nGreenValue,
569 : nBlueValue, nDstNoData));
570 19399 : iIndex = *psIndex;
571 : }
572 : else
573 : {
574 60601 : iIndex = *psIndex;
575 : }
576 : }
577 :
578 1474960 : pabyIndex[i] = static_cast<GByte>(iIndex);
579 1474960 : if (!bDither)
580 287500 : continue;
581 :
582 : /* --------------------------------------------------------------------
583 : */
584 : /* Compute Red error, and carry it on to the next error line.
585 : */
586 : /* --------------------------------------------------------------------
587 : */
588 1187460 : nError = nRedValue - CAST_PCT(anPCT)[4 * iIndex + 0];
589 1187460 : nSixth = nError / 6;
590 :
591 1187460 : panError[i * 3] += nSixth;
592 1187460 : panError[i * 3 + 6] = nSixth;
593 1187460 : panError[i * 3 + 3] += nError - 5 * nSixth;
594 :
595 1187460 : nLastRedError = 2 * nSixth;
596 :
597 : /* --------------------------------------------------------------------
598 : */
599 : /* Compute Green error, and carry it on to the next error line.
600 : */
601 : /* --------------------------------------------------------------------
602 : */
603 1187460 : nError = nGreenValue - CAST_PCT(anPCT)[4 * iIndex + 1];
604 1187460 : nSixth = nError / 6;
605 :
606 1187460 : panError[i * 3 + 1] += nSixth;
607 1187460 : panError[i * 3 + 6 + 1] = nSixth;
608 1187460 : panError[i * 3 + 3 + 1] += nError - 5 * nSixth;
609 :
610 1187460 : nLastGreenError = 2 * nSixth;
611 :
612 : /* --------------------------------------------------------------------
613 : */
614 : /* Compute Blue error, and carry it on to the next error line.
615 : */
616 : /* --------------------------------------------------------------------
617 : */
618 1187460 : nError = nBlueValue - CAST_PCT(anPCT)[4 * iIndex + 2];
619 1187460 : nSixth = nError / 6;
620 :
621 1187460 : panError[i * 3 + 2] += nSixth;
622 1187460 : panError[i * 3 + 6 + 2] = nSixth;
623 1187460 : panError[i * 3 + 3 + 2] += nError - 5 * nSixth;
624 :
625 1187460 : nLastBlueError = 2 * nSixth;
626 : }
627 :
628 : /* --------------------------------------------------------------------
629 : */
630 : /* Write results. */
631 : /* --------------------------------------------------------------------
632 : */
633 5174 : err = GDALRasterIO(hTarget, GF_Write, 0, iScanline, nXSize, 1,
634 : pabyIndex, nXSize, 1, GDT_UInt8, 0, 0);
635 5174 : if (err != CE_None)
636 0 : break;
637 : }
638 :
639 32 : pfnProgress(1.0, nullptr, pProgressArg);
640 :
641 : /* -------------------------------------------------------------------- */
642 : /* Cleanup */
643 : /* -------------------------------------------------------------------- */
644 32 : CPLFree(pabyRed);
645 32 : CPLFree(pabyGreen);
646 32 : CPLFree(pabyBlue);
647 32 : CPLFree(pabyIndex);
648 32 : CPLFree(panError);
649 32 : CPLFree(pabyColorMap);
650 :
651 32 : return err;
652 : }
653 :
654 875508 : static int FindNearestColor(int nColors, const int *panPCT, int nRedValue,
655 : int nGreenValue, int nBlueValue, int nDstNoData)
656 :
657 : {
658 : #ifdef USE_SSE2
659 875508 : int nBestDist = 768;
660 875508 : int nBestIndex = 0;
661 :
662 875508 : int anDistanceUnaligned[16 + 4] =
663 : {}; // 4 for alignment on 16-byte boundary.
664 875508 : int *anDistance = ALIGN_INT_ARRAY_ON_16_BYTE(anDistanceUnaligned);
665 :
666 875508 : const __m128i ff = _mm_set1_epi32(0xFFFFFFFF);
667 875508 : const __m128i mask_low = _mm_srli_epi64(ff, 32);
668 875508 : const __m128i mask_high = _mm_slli_epi64(ff, 32);
669 :
670 : const unsigned int nColorVal =
671 875508 : MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue);
672 1751020 : const __m128i thisColor = _mm_set1_epi32(nColorVal);
673 875508 : const __m128i thisColor_low = _mm_srli_epi64(thisColor, 32);
674 875508 : const __m128i thisColor_high = _mm_slli_epi64(thisColor, 32);
675 :
676 26892900 : for (int iColor = 0; iColor < nColors; iColor += 8)
677 : {
678 : const __m128i pctColor =
679 26017400 : _mm_load_si128(reinterpret_cast<const __m128i *>(&panPCT[iColor]));
680 26017400 : const __m128i pctColor2 = _mm_load_si128(
681 26017400 : reinterpret_cast<const __m128i *>(&panPCT[iColor + 4]));
682 :
683 52034800 : _mm_store_si128(
684 : reinterpret_cast<__m128i *>(anDistance),
685 : _mm_sad_epu8(_mm_and_si128(pctColor, mask_low), thisColor_low));
686 26017400 : _mm_store_si128(
687 26017400 : reinterpret_cast<__m128i *>(anDistance + 4),
688 : _mm_sad_epu8(_mm_and_si128(pctColor, mask_high), thisColor_high));
689 26017400 : _mm_store_si128(
690 26017400 : reinterpret_cast<__m128i *>(anDistance + 8),
691 : _mm_sad_epu8(_mm_and_si128(pctColor2, mask_low), thisColor_low));
692 26017400 : _mm_store_si128(
693 26017400 : reinterpret_cast<__m128i *>(anDistance + 12),
694 : _mm_sad_epu8(_mm_and_si128(pctColor2, mask_high), thisColor_high));
695 :
696 26017400 : if (anDistance[0] < nBestDist && iColor != nDstNoData)
697 : {
698 1119370 : nBestIndex = iColor;
699 1119370 : nBestDist = anDistance[0];
700 : }
701 26017400 : if (anDistance[4] < nBestDist && iColor + 1 != nDstNoData)
702 : {
703 800944 : nBestIndex = iColor + 1;
704 800944 : nBestDist = anDistance[4];
705 : }
706 26017400 : if (anDistance[2] < nBestDist && iColor + 2 != nDstNoData)
707 : {
708 517890 : nBestIndex = iColor + 2;
709 517890 : nBestDist = anDistance[2];
710 : }
711 26017400 : if (anDistance[6] < nBestDist && iColor + 3 != nDstNoData)
712 : {
713 494934 : nBestIndex = iColor + 3;
714 494934 : nBestDist = anDistance[6];
715 : }
716 26017400 : if (anDistance[8 + 0] < nBestDist && iColor + 4 != nDstNoData)
717 : {
718 363062 : nBestIndex = iColor + 4;
719 363062 : nBestDist = anDistance[8 + 0];
720 : }
721 26017400 : if (anDistance[8 + 4] < nBestDist && iColor + 5 != nDstNoData)
722 : {
723 318223 : nBestIndex = iColor + 4 + 1;
724 318223 : nBestDist = anDistance[8 + 4];
725 : }
726 26017400 : if (anDistance[8 + 2] < nBestDist && iColor + 6 != nDstNoData)
727 : {
728 312654 : nBestIndex = iColor + 4 + 2;
729 312654 : nBestDist = anDistance[8 + 2];
730 : }
731 26017400 : if (anDistance[8 + 6] < nBestDist && iColor + 7 != nDstNoData)
732 : {
733 393373 : nBestIndex = iColor + 4 + 3;
734 393373 : nBestDist = anDistance[8 + 6];
735 : }
736 : }
737 875508 : return nBestIndex;
738 : #else
739 : int nBestDist = 768;
740 : int nBestIndex = 0;
741 :
742 : for (int iColor = 0; iColor < nColors; iColor++)
743 : {
744 : if (iColor != nDstNoData)
745 : {
746 : const int nThisDist =
747 : std::abs(nRedValue - panPCT[4 * iColor + 0]) +
748 : std::abs(nGreenValue - panPCT[4 * iColor + 1]) +
749 : std::abs(nBlueValue - panPCT[4 * iColor + 2]);
750 :
751 : if (nThisDist < nBestDist)
752 : {
753 : nBestIndex = iColor;
754 : nBestDist = nThisDist;
755 : }
756 : }
757 : }
758 : return nBestIndex;
759 : #endif
760 : }
761 :
762 : /************************************************************************/
763 : /* FindNearestColor() */
764 : /* */
765 : /* Finear near PCT color for any RGB color. */
766 : /************************************************************************/
767 :
768 22 : static void FindNearestColor(int nColors, const int *panPCT,
769 : GByte *pabyColorMap, int nCLevels, int nDstNoData)
770 :
771 : {
772 : /* -------------------------------------------------------------------- */
773 : /* Loop over all the cells in the high density cube. */
774 : /* -------------------------------------------------------------------- */
775 726 : for (int iBlue = 0; iBlue < nCLevels; iBlue++)
776 : {
777 23232 : for (int iGreen = 0; iGreen < nCLevels; iGreen++)
778 : {
779 743424 : for (int iRed = 0; iRed < nCLevels; iRed++)
780 : {
781 720896 : const int nRedValue = (iRed * 255) / (nCLevels - 1);
782 720896 : const int nGreenValue = (iGreen * 255) / (nCLevels - 1);
783 720896 : const int nBlueValue = (iBlue * 255) / (nCLevels - 1);
784 :
785 : const int nBestIndex =
786 720896 : FindNearestColor(nColors, panPCT, nRedValue, nGreenValue,
787 : nBlueValue, nDstNoData);
788 720896 : pabyColorMap[iRed + iGreen * nCLevels +
789 720896 : iBlue * nCLevels * nCLevels] =
790 : static_cast<GByte>(nBestIndex);
791 : }
792 : }
793 : }
794 22 : }
|
