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