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 <cstdlib>
38 : #include <cstring>
39 : #include <algorithm>
40 :
41 : #include "cpl_conv.h"
42 : #include "cpl_error.h"
43 : #include "cpl_progress.h"
44 : #include "cpl_vsi.h"
45 : #include "gdal.h"
46 : #include "gdal_priv.h"
47 :
48 : #ifdef USE_NEON_OPTIMIZATIONS
49 : #define USE_SSE2
50 : #include "include_sse2neon.h"
51 : #elif defined(__x86_64) || defined(_M_X64)
52 : #define USE_SSE2
53 : #include <emmintrin.h>
54 : #endif
55 :
56 : #ifdef USE_SSE2
57 :
58 : #define CAST_PCT(x) reinterpret_cast<GByte *>(x)
59 : #define ALIGN_INT_ARRAY_ON_16_BYTE(x) \
60 : (((reinterpret_cast<GUIntptr_t>(x) % 16) != 0) \
61 : ? reinterpret_cast<int *>(reinterpret_cast<GByte *>(x) + 16 - \
62 : (reinterpret_cast<GUIntptr_t>(x) % 16)) \
63 : : (x))
64 : #else
65 : #define CAST_PCT(x) x
66 : #endif
67 :
68 : #ifndef MAKE_COLOR_CODE_defined
69 : #define MAKE_COLOR_CODE_defined
70 :
71 820864 : static int MAKE_COLOR_CODE(int r, int g, int b)
72 : {
73 820864 : return r | (g << 8) | (b << 16);
74 : }
75 : #endif
76 :
77 : static void FindNearestColor(int nColors, int *panPCT, GByte *pabyColorMap,
78 : int nCLevels);
79 : static int FindNearestColor(int nColors, int *panPCT, int nRedValue,
80 : int nGreenValue, int nBlueValue);
81 :
82 : // Structure for a hashmap from a color code to a color index of the
83 : // color table.
84 :
85 : // NOTE: if changing the size of this structure, edit
86 : // MEDIAN_CUT_AND_DITHER_BUFFER_SIZE_65536 in gdal_alg_priv.h and take
87 : // into account HashHistogram in gdalmediancut.cpp.
88 : typedef struct
89 : {
90 : GUInt32 nColorCode;
91 : GUInt32 nColorCode2;
92 : GUInt32 nColorCode3;
93 : GByte nIndex;
94 : GByte nIndex2;
95 : GByte nIndex3;
96 : GByte nPadding;
97 : } ColorIndex;
98 :
99 : /************************************************************************/
100 : /* GDALDitherRGB2PCT() */
101 : /************************************************************************/
102 :
103 : #ifndef IsColorCodeSet_defined
104 : #define IsColorCodeSet_defined
105 :
106 146386 : static inline bool IsColorCodeSet(GUInt32 nColorCode)
107 : {
108 146386 : return (nColorCode >> 31) == 0;
109 : }
110 : #endif
111 :
112 : /**
113 : * 24bit to 8bit conversion with dithering.
114 : *
115 : * This functions utilizes Floyd-Steinberg dithering in the process of
116 : * converting a 24bit RGB image into a pseudocolored 8bit image using a
117 : * provided color table.
118 : *
119 : * The red, green and blue input bands do not necessarily need to come
120 : * from the same file, but they must be the same width and height. They will
121 : * be clipped to 8bit during reading, so non-eight bit bands are generally
122 : * inappropriate. Likewise the hTarget band will be written with 8bit values
123 : * and must match the width and height of the source bands.
124 : *
125 : * The color table cannot have more than 256 entries.
126 : *
127 : * @param hRed Red input band.
128 : * @param hGreen Green input band.
129 : * @param hBlue Blue input band.
130 : * @param hTarget Output band.
131 : * @param hColorTable the color table to use with the output band.
132 : * @param pfnProgress callback for reporting algorithm progress matching the
133 : * GDALProgressFunc() semantics. May be NULL.
134 : * @param pProgressArg callback argument passed to pfnProgress.
135 : *
136 : * @return CE_None on success or CE_Failure if an error occurs.
137 : */
138 :
139 6 : int CPL_STDCALL GDALDitherRGB2PCT(GDALRasterBandH hRed, GDALRasterBandH hGreen,
140 : GDALRasterBandH hBlue,
141 : GDALRasterBandH hTarget,
142 : GDALColorTableH hColorTable,
143 : GDALProgressFunc pfnProgress,
144 : void *pProgressArg)
145 :
146 : {
147 6 : return GDALDitherRGB2PCTInternal(hRed, hGreen, hBlue, hTarget, hColorTable,
148 : 5, nullptr, TRUE, pfnProgress,
149 6 : pProgressArg);
150 : }
151 :
152 16 : int GDALDitherRGB2PCTInternal(
153 : GDALRasterBandH hRed, GDALRasterBandH hGreen, GDALRasterBandH hBlue,
154 : GDALRasterBandH hTarget, GDALColorTableH hColorTable, int nBits,
155 : // NULL or at least 256 * 256 * 256 * sizeof(GInt16) bytes.
156 : GInt16 *pasDynamicColorMap, int bDither, GDALProgressFunc pfnProgress,
157 : void *pProgressArg)
158 : {
159 16 : VALIDATE_POINTER1(hRed, "GDALDitherRGB2PCT", CE_Failure);
160 16 : VALIDATE_POINTER1(hGreen, "GDALDitherRGB2PCT", CE_Failure);
161 16 : VALIDATE_POINTER1(hBlue, "GDALDitherRGB2PCT", CE_Failure);
162 16 : VALIDATE_POINTER1(hTarget, "GDALDitherRGB2PCT", CE_Failure);
163 16 : VALIDATE_POINTER1(hColorTable, "GDALDitherRGB2PCT", CE_Failure);
164 :
165 : /* -------------------------------------------------------------------- */
166 : /* Validate parameters. */
167 : /* -------------------------------------------------------------------- */
168 16 : const int nXSize = GDALGetRasterBandXSize(hRed);
169 16 : const int nYSize = GDALGetRasterBandYSize(hRed);
170 :
171 16 : if (GDALGetRasterBandXSize(hGreen) != nXSize ||
172 16 : GDALGetRasterBandYSize(hGreen) != nYSize ||
173 48 : GDALGetRasterBandXSize(hBlue) != nXSize ||
174 16 : GDALGetRasterBandYSize(hBlue) != nYSize)
175 : {
176 0 : CPLError(CE_Failure, CPLE_IllegalArg,
177 : "Green or blue band doesn't match size of red band.");
178 :
179 0 : return CE_Failure;
180 : }
181 :
182 32 : if (GDALGetRasterBandXSize(hTarget) != nXSize ||
183 16 : GDALGetRasterBandYSize(hTarget) != nYSize)
184 : {
185 0 : CPLError(CE_Failure, CPLE_IllegalArg,
186 : "GDALDitherRGB2PCT(): "
187 : "Target band doesn't match size of source bands.");
188 :
189 0 : return CE_Failure;
190 : }
191 :
192 16 : if (pfnProgress == nullptr)
193 11 : pfnProgress = GDALDummyProgress;
194 :
195 : /* -------------------------------------------------------------------- */
196 : /* Setup more direct colormap. */
197 : /* -------------------------------------------------------------------- */
198 : int iColor;
199 : #ifdef USE_SSE2
200 : int anPCTUnaligned[256 + 4]; // 4 for alignment on 16-byte boundary.
201 16 : int *anPCT = ALIGN_INT_ARRAY_ON_16_BYTE(anPCTUnaligned);
202 : #else
203 : int anPCT[256 * 4] = {};
204 : #endif
205 16 : const int nColors = GDALGetColorEntryCount(hColorTable);
206 :
207 16 : if (nColors == 0)
208 : {
209 0 : CPLError(CE_Failure, CPLE_IllegalArg,
210 : "GDALDitherRGB2PCT(): "
211 : "Color table must not be empty.");
212 :
213 0 : return CE_Failure;
214 : }
215 16 : else if (nColors > 256)
216 : {
217 0 : CPLError(CE_Failure, CPLE_IllegalArg,
218 : "GDALDitherRGB2PCT(): "
219 : "Color table cannot have more than 256 entries.");
220 :
221 0 : return CE_Failure;
222 : }
223 :
224 16 : iColor = 0;
225 3337 : do
226 : {
227 : GDALColorEntry sEntry;
228 :
229 3353 : GDALGetColorEntryAsRGB(hColorTable, iColor, &sEntry);
230 3353 : CAST_PCT(anPCT)[4 * iColor + 0] = static_cast<GByte>(sEntry.c1);
231 3353 : CAST_PCT(anPCT)[4 * iColor + 1] = static_cast<GByte>(sEntry.c2);
232 3353 : CAST_PCT(anPCT)[4 * iColor + 2] = static_cast<GByte>(sEntry.c3);
233 3353 : CAST_PCT(anPCT)[4 * iColor + 3] = 0;
234 :
235 3353 : iColor++;
236 3353 : } while (iColor < nColors);
237 :
238 : #ifdef USE_SSE2
239 : // Pad to multiple of 8 colors.
240 16 : const int nColorsMod8 = nColors % 8;
241 16 : if (nColorsMod8)
242 : {
243 1 : int iDest = nColors;
244 8 : for (iColor = 0; iColor < 8 - nColorsMod8 && iDest < 256;
245 : iColor++, iDest++)
246 : {
247 7 : anPCT[iDest] = anPCT[nColors - 1];
248 : }
249 : }
250 : #endif
251 :
252 : /* -------------------------------------------------------------------- */
253 : /* Setup various variables. */
254 : /* -------------------------------------------------------------------- */
255 16 : const int nCLevels = 1 << nBits;
256 16 : const int nCLevelsCube = nCLevels * nCLevels * nCLevels;
257 16 : ColorIndex *psColorIndexMap = nullptr;
258 :
259 16 : GByte *pabyRed = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
260 16 : GByte *pabyGreen = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
261 16 : GByte *pabyBlue = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
262 :
263 16 : GByte *pabyIndex = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize));
264 :
265 : int *panError =
266 16 : static_cast<int *>(VSI_CALLOC_VERBOSE(sizeof(int), (nXSize + 2) * 3));
267 :
268 16 : if (pabyRed == nullptr || pabyGreen == nullptr || pabyBlue == nullptr ||
269 16 : pabyIndex == nullptr || panError == nullptr)
270 : {
271 0 : CPLFree(pabyRed);
272 0 : CPLFree(pabyGreen);
273 0 : CPLFree(pabyBlue);
274 0 : CPLFree(pabyIndex);
275 0 : CPLFree(panError);
276 :
277 0 : return CE_Failure;
278 : }
279 :
280 16 : GByte *pabyColorMap = nullptr;
281 16 : if (pasDynamicColorMap == nullptr)
282 : {
283 : /* --------------------------------------------------------------------
284 : */
285 : /* Build a 24bit to 8 bit color mapping. */
286 : /* --------------------------------------------------------------------
287 : */
288 :
289 : pabyColorMap = static_cast<GByte *>(
290 6 : VSI_MALLOC_VERBOSE(nCLevelsCube * sizeof(GByte)));
291 6 : if (pabyColorMap == nullptr)
292 : {
293 0 : CPLFree(pabyRed);
294 0 : CPLFree(pabyGreen);
295 0 : CPLFree(pabyBlue);
296 0 : CPLFree(pabyIndex);
297 0 : CPLFree(panError);
298 0 : CPLFree(pabyColorMap);
299 :
300 0 : return CE_Failure;
301 : }
302 :
303 6 : FindNearestColor(nColors, anPCT, pabyColorMap, nCLevels);
304 : }
305 : else
306 : {
307 10 : pabyColorMap = nullptr;
308 10 : if (nBits == 8 && static_cast<GIntBig>(nXSize) * nYSize <= 65536)
309 : {
310 : // If the image is small enough, then the number of colors
311 : // will be limited and using a hashmap, rather than a full table
312 : // will be more efficient.
313 9 : psColorIndexMap =
314 : reinterpret_cast<ColorIndex *>(pasDynamicColorMap);
315 9 : memset(psColorIndexMap, 0xFF, sizeof(ColorIndex) * PRIME_FOR_65536);
316 : }
317 : else
318 : {
319 1 : memset(pasDynamicColorMap, 0xFF, 256 * 256 * 256 * sizeof(GInt16));
320 : }
321 : }
322 :
323 : /* ==================================================================== */
324 : /* Loop over all scanlines of data to process. */
325 : /* ==================================================================== */
326 16 : CPLErr err = CE_None;
327 :
328 2290 : for (int iScanline = 0; iScanline < nYSize; iScanline++)
329 : {
330 : /* --------------------------------------------------------------------
331 : */
332 : /* Report progress */
333 : /* --------------------------------------------------------------------
334 : */
335 2274 : if (!pfnProgress(iScanline / static_cast<double>(nYSize), nullptr,
336 : pProgressArg))
337 : {
338 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User Terminated");
339 0 : CPLFree(pabyRed);
340 0 : CPLFree(pabyGreen);
341 0 : CPLFree(pabyBlue);
342 0 : CPLFree(pabyIndex);
343 0 : CPLFree(panError);
344 0 : CPLFree(pabyColorMap);
345 :
346 0 : return CE_Failure;
347 : }
348 :
349 : /* --------------------------------------------------------------------
350 : */
351 : /* Read source data. */
352 : /* --------------------------------------------------------------------
353 : */
354 2274 : CPLErr err1 = GDALRasterIO(hRed, GF_Read, 0, iScanline, nXSize, 1,
355 : pabyRed, nXSize, 1, GDT_Byte, 0, 0);
356 2274 : if (err1 == CE_None)
357 2274 : err1 = GDALRasterIO(hGreen, GF_Read, 0, iScanline, nXSize, 1,
358 : pabyGreen, nXSize, 1, GDT_Byte, 0, 0);
359 2274 : if (err1 == CE_None)
360 2274 : err1 = GDALRasterIO(hBlue, GF_Read, 0, iScanline, nXSize, 1,
361 : pabyBlue, nXSize, 1, GDT_Byte, 0, 0);
362 2274 : if (err1 != CE_None)
363 : {
364 0 : CPLFree(pabyRed);
365 0 : CPLFree(pabyGreen);
366 0 : CPLFree(pabyBlue);
367 0 : CPLFree(pabyIndex);
368 0 : CPLFree(panError);
369 0 : CPLFree(pabyColorMap);
370 :
371 0 : return err1;
372 : }
373 :
374 : /* --------------------------------------------------------------------
375 : */
376 : /* Apply the error from the previous line to this one. */
377 : /* --------------------------------------------------------------------
378 : */
379 2274 : if (bDither)
380 : {
381 278468 : for (int i = 0; i < nXSize; i++)
382 : {
383 277144 : pabyRed[i] = static_cast<GByte>(std::max(
384 277144 : 0, std::min(255, (pabyRed[i] + panError[i * 3 + 0 + 3]))));
385 277144 : pabyGreen[i] = static_cast<GByte>(std::max(
386 554288 : 0,
387 277144 : std::min(255, (pabyGreen[i] + panError[i * 3 + 1 + 3]))));
388 277144 : pabyBlue[i] = static_cast<GByte>(std::max(
389 277144 : 0, std::min(255, (pabyBlue[i] + panError[i * 3 + 2 + 3]))));
390 : }
391 :
392 1324 : memset(panError, 0, sizeof(int) * (nXSize + 2) * 3);
393 : }
394 :
395 : /* --------------------------------------------------------------------
396 : */
397 : /* Figure out the nearest color to the RGB value. */
398 : /* --------------------------------------------------------------------
399 : */
400 2274 : int nLastRedError = 0;
401 2274 : int nLastGreenError = 0;
402 2274 : int nLastBlueError = 0;
403 :
404 486918 : for (int i = 0; i < nXSize; i++)
405 : {
406 : const int nRedValue =
407 484644 : std::max(0, std::min(255, pabyRed[i] + nLastRedError));
408 : const int nGreenValue =
409 484644 : std::max(0, std::min(255, pabyGreen[i] + nLastGreenError));
410 : const int nBlueValue =
411 484644 : std::max(0, std::min(255, pabyBlue[i] + nLastBlueError));
412 :
413 484644 : int iIndex = 0;
414 484644 : int nError = 0;
415 484644 : int nSixth = 0;
416 484644 : if (psColorIndexMap)
417 : {
418 : const GUInt32 nColorCode =
419 389644 : MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue);
420 389644 : GUInt32 nIdx = nColorCode % PRIME_FOR_65536;
421 : while (true)
422 : {
423 389651 : if (psColorIndexMap[nIdx].nColorCode == nColorCode)
424 : {
425 252924 : iIndex = psColorIndexMap[nIdx].nIndex;
426 252924 : break;
427 : }
428 136727 : if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode))
429 : {
430 125973 : psColorIndexMap[nIdx].nColorCode = nColorCode;
431 125973 : iIndex = FindNearestColor(nColors, anPCT, nRedValue,
432 : nGreenValue, nBlueValue);
433 125973 : psColorIndexMap[nIdx].nIndex =
434 : static_cast<GByte>(iIndex);
435 125973 : break;
436 : }
437 10754 : if (psColorIndexMap[nIdx].nColorCode2 == nColorCode)
438 : {
439 1476 : iIndex = psColorIndexMap[nIdx].nIndex2;
440 1476 : break;
441 : }
442 9278 : if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode2))
443 : {
444 8876 : psColorIndexMap[nIdx].nColorCode2 = nColorCode;
445 8876 : iIndex = FindNearestColor(nColors, anPCT, nRedValue,
446 : nGreenValue, nBlueValue);
447 8876 : psColorIndexMap[nIdx].nIndex2 =
448 : static_cast<GByte>(iIndex);
449 8876 : break;
450 : }
451 402 : if (psColorIndexMap[nIdx].nColorCode3 == nColorCode)
452 : {
453 31 : iIndex = psColorIndexMap[nIdx].nIndex3;
454 31 : break;
455 : }
456 371 : if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode3))
457 : {
458 364 : psColorIndexMap[nIdx].nColorCode3 = nColorCode;
459 364 : iIndex = FindNearestColor(nColors, anPCT, nRedValue,
460 : nGreenValue, nBlueValue);
461 364 : psColorIndexMap[nIdx].nIndex3 =
462 : static_cast<GByte>(iIndex);
463 364 : break;
464 : }
465 :
466 0 : do
467 : {
468 7 : nIdx += 257;
469 7 : if (nIdx >= PRIME_FOR_65536)
470 0 : nIdx -= PRIME_FOR_65536;
471 : } while (
472 7 : IsColorCodeSet(psColorIndexMap[nIdx].nColorCode) &&
473 6 : psColorIndexMap[nIdx].nColorCode != nColorCode &&
474 3 : IsColorCodeSet(psColorIndexMap[nIdx].nColorCode2) &&
475 0 : psColorIndexMap[nIdx].nColorCode2 != nColorCode &&
476 10 : IsColorCodeSet(psColorIndexMap[nIdx].nColorCode3) &&
477 0 : psColorIndexMap[nIdx].nColorCode3 != nColorCode);
478 : }
479 : }
480 95000 : else if (pasDynamicColorMap == nullptr)
481 : {
482 15000 : const int iRed = nRedValue * nCLevels / 256;
483 15000 : const int iGreen = nGreenValue * nCLevels / 256;
484 15000 : const int iBlue = nBlueValue * nCLevels / 256;
485 :
486 15000 : iIndex = pabyColorMap[iRed + iGreen * nCLevels +
487 15000 : iBlue * nCLevels * nCLevels];
488 : }
489 : else
490 : {
491 : const GUInt32 nColorCode =
492 80000 : MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue);
493 80000 : GInt16 *psIndex = &pasDynamicColorMap[nColorCode];
494 80000 : if (*psIndex < 0)
495 : {
496 19399 : *psIndex = static_cast<GInt16>(FindNearestColor(
497 : nColors, anPCT, nRedValue, nGreenValue, nBlueValue));
498 19399 : iIndex = *psIndex;
499 : }
500 : else
501 : {
502 60601 : iIndex = *psIndex;
503 : }
504 : }
505 :
506 484644 : pabyIndex[i] = static_cast<GByte>(iIndex);
507 484644 : if (!bDither)
508 207500 : continue;
509 :
510 : /* --------------------------------------------------------------------
511 : */
512 : /* Compute Red error, and carry it on to the next error line.
513 : */
514 : /* --------------------------------------------------------------------
515 : */
516 277144 : nError = nRedValue - CAST_PCT(anPCT)[4 * iIndex + 0];
517 277144 : nSixth = nError / 6;
518 :
519 277144 : panError[i * 3] += nSixth;
520 277144 : panError[i * 3 + 6] = nSixth;
521 277144 : panError[i * 3 + 3] += nError - 5 * nSixth;
522 :
523 277144 : nLastRedError = 2 * nSixth;
524 :
525 : /* --------------------------------------------------------------------
526 : */
527 : /* Compute Green error, and carry it on to the next error line.
528 : */
529 : /* --------------------------------------------------------------------
530 : */
531 277144 : nError = nGreenValue - CAST_PCT(anPCT)[4 * iIndex + 1];
532 277144 : nSixth = nError / 6;
533 :
534 277144 : panError[i * 3 + 1] += nSixth;
535 277144 : panError[i * 3 + 6 + 1] = nSixth;
536 277144 : panError[i * 3 + 3 + 1] += nError - 5 * nSixth;
537 :
538 277144 : nLastGreenError = 2 * nSixth;
539 :
540 : /* --------------------------------------------------------------------
541 : */
542 : /* Compute Blue error, and carry it on to the next error line.
543 : */
544 : /* --------------------------------------------------------------------
545 : */
546 277144 : nError = nBlueValue - CAST_PCT(anPCT)[4 * iIndex + 2];
547 277144 : nSixth = nError / 6;
548 :
549 277144 : panError[i * 3 + 2] += nSixth;
550 277144 : panError[i * 3 + 6 + 2] = nSixth;
551 277144 : panError[i * 3 + 3 + 2] += nError - 5 * nSixth;
552 :
553 277144 : nLastBlueError = 2 * nSixth;
554 : }
555 :
556 : /* --------------------------------------------------------------------
557 : */
558 : /* Write results. */
559 : /* --------------------------------------------------------------------
560 : */
561 2274 : err = GDALRasterIO(hTarget, GF_Write, 0, iScanline, nXSize, 1,
562 : pabyIndex, nXSize, 1, GDT_Byte, 0, 0);
563 2274 : if (err != CE_None)
564 0 : break;
565 : }
566 :
567 16 : pfnProgress(1.0, nullptr, pProgressArg);
568 :
569 : /* -------------------------------------------------------------------- */
570 : /* Cleanup */
571 : /* -------------------------------------------------------------------- */
572 16 : CPLFree(pabyRed);
573 16 : CPLFree(pabyGreen);
574 16 : CPLFree(pabyBlue);
575 16 : CPLFree(pabyIndex);
576 16 : CPLFree(panError);
577 16 : CPLFree(pabyColorMap);
578 :
579 16 : return err;
580 : }
581 :
582 351220 : static int FindNearestColor(int nColors, int *panPCT, int nRedValue,
583 : int nGreenValue, int nBlueValue)
584 :
585 : {
586 : #ifdef USE_SSE2
587 351220 : int nBestDist = 768;
588 351220 : int nBestIndex = 0;
589 :
590 351220 : int anDistanceUnaligned[16 + 4] =
591 : {}; // 4 for alignment on 16-byte boundary.
592 351220 : int *anDistance = ALIGN_INT_ARRAY_ON_16_BYTE(anDistanceUnaligned);
593 :
594 351220 : const __m128i ff = _mm_set1_epi32(0xFFFFFFFF);
595 351220 : const __m128i mask_low = _mm_srli_epi64(ff, 32);
596 351220 : const __m128i mask_high = _mm_slli_epi64(ff, 32);
597 :
598 : const unsigned int nColorVal =
599 351220 : MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue);
600 702440 : const __m128i thisColor = _mm_set1_epi32(nColorVal);
601 351220 : const __m128i thisColor_low = _mm_srli_epi64(thisColor, 32);
602 351220 : const __m128i thisColor_high = _mm_slli_epi64(thisColor, 32);
603 :
604 9591380 : for (int iColor = 0; iColor < nColors; iColor += 8)
605 : {
606 : const __m128i pctColor =
607 9240160 : _mm_load_si128(reinterpret_cast<__m128i *>(&panPCT[iColor]));
608 : const __m128i pctColor2 =
609 18480300 : _mm_load_si128(reinterpret_cast<__m128i *>(&panPCT[iColor + 4]));
610 :
611 18480300 : _mm_store_si128(
612 : reinterpret_cast<__m128i *>(anDistance),
613 : _mm_sad_epu8(_mm_and_si128(pctColor, mask_low), thisColor_low));
614 9240160 : _mm_store_si128(
615 9240160 : reinterpret_cast<__m128i *>(anDistance + 4),
616 : _mm_sad_epu8(_mm_and_si128(pctColor, mask_high), thisColor_high));
617 9240160 : _mm_store_si128(
618 9240160 : reinterpret_cast<__m128i *>(anDistance + 8),
619 : _mm_sad_epu8(_mm_and_si128(pctColor2, mask_low), thisColor_low));
620 9240160 : _mm_store_si128(
621 9240160 : reinterpret_cast<__m128i *>(anDistance + 12),
622 : _mm_sad_epu8(_mm_and_si128(pctColor2, mask_high), thisColor_high));
623 :
624 9240160 : if (anDistance[0] < nBestDist)
625 : {
626 441401 : nBestIndex = iColor;
627 441401 : nBestDist = anDistance[0];
628 : }
629 9240160 : if (anDistance[4] < nBestDist)
630 : {
631 271510 : nBestIndex = iColor + 1;
632 271510 : nBestDist = anDistance[4];
633 : }
634 9240160 : if (anDistance[2] < nBestDist)
635 : {
636 150984 : nBestIndex = iColor + 2;
637 150984 : nBestDist = anDistance[2];
638 : }
639 9240160 : if (anDistance[6] < nBestDist)
640 : {
641 175743 : nBestIndex = iColor + 3;
642 175743 : nBestDist = anDistance[6];
643 : }
644 9240160 : if (anDistance[8 + 0] < nBestDist)
645 : {
646 119903 : nBestIndex = iColor + 4;
647 119903 : nBestDist = anDistance[8 + 0];
648 : }
649 9240160 : if (anDistance[8 + 4] < nBestDist)
650 : {
651 98139 : nBestIndex = iColor + 4 + 1;
652 98139 : nBestDist = anDistance[8 + 4];
653 : }
654 9240160 : if (anDistance[8 + 2] < nBestDist)
655 : {
656 132829 : nBestIndex = iColor + 4 + 2;
657 132829 : nBestDist = anDistance[8 + 2];
658 : }
659 9240160 : if (anDistance[8 + 6] < nBestDist)
660 : {
661 185338 : nBestIndex = iColor + 4 + 3;
662 185338 : nBestDist = anDistance[8 + 6];
663 : }
664 : }
665 351220 : return nBestIndex;
666 : #else
667 : int nBestDist = 768;
668 : int nBestIndex = 0;
669 :
670 : for (int iColor = 0; iColor < nColors; iColor++)
671 : {
672 : const int nThisDist = std::abs(nRedValue - panPCT[4 * iColor + 0]) +
673 : std::abs(nGreenValue - panPCT[4 * iColor + 1]) +
674 : std::abs(nBlueValue - panPCT[4 * iColor + 2]);
675 :
676 : if (nThisDist < nBestDist)
677 : {
678 : nBestIndex = iColor;
679 : nBestDist = nThisDist;
680 : }
681 : }
682 : return nBestIndex;
683 : #endif
684 : }
685 :
686 : /************************************************************************/
687 : /* FindNearestColor() */
688 : /* */
689 : /* Finear near PCT color for any RGB color. */
690 : /************************************************************************/
691 :
692 6 : static void FindNearestColor(int nColors, int *panPCT, GByte *pabyColorMap,
693 : int nCLevels)
694 :
695 : {
696 : /* -------------------------------------------------------------------- */
697 : /* Loop over all the cells in the high density cube. */
698 : /* -------------------------------------------------------------------- */
699 198 : for (int iBlue = 0; iBlue < nCLevels; iBlue++)
700 : {
701 6336 : for (int iGreen = 0; iGreen < nCLevels; iGreen++)
702 : {
703 202752 : for (int iRed = 0; iRed < nCLevels; iRed++)
704 : {
705 196608 : const int nRedValue = (iRed * 255) / (nCLevels - 1);
706 196608 : const int nGreenValue = (iGreen * 255) / (nCLevels - 1);
707 196608 : const int nBlueValue = (iBlue * 255) / (nCLevels - 1);
708 :
709 196608 : const int nBestIndex = FindNearestColor(
710 : nColors, panPCT, nRedValue, nGreenValue, nBlueValue);
711 196608 : pabyColorMap[iRed + iGreen * nCLevels +
712 196608 : iBlue * nCLevels * nCLevels] =
713 : static_cast<GByte>(nBestIndex);
714 : }
715 : }
716 : }
717 6 : }
|
