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