/src/gdal/alg/gdaldither.cpp
Line | Count | Source (jump to first uncovered line) |
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 | 0 | #define CAST_PCT(x) reinterpret_cast<GByte *>(x) |
59 | | #define ALIGN_INT_ARRAY_ON_16_BYTE(x) \ |
60 | 0 | (((reinterpret_cast<GUIntptr_t>(x) % 16) != 0) \ |
61 | 0 | ? reinterpret_cast<int *>(reinterpret_cast<GByte *>(x) + 16 - \ |
62 | 0 | (reinterpret_cast<GUIntptr_t>(x) % 16)) \ |
63 | 0 | : (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 | | static int MAKE_COLOR_CODE(int r, int g, int b) |
72 | 0 | { |
73 | 0 | return r | (g << 8) | (b << 16); |
74 | 0 | } |
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 | | static inline bool IsColorCodeSet(GUInt32 nColorCode) |
107 | 0 | { |
108 | 0 | return (nColorCode >> 31) == 0; |
109 | 0 | } |
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 | | 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 | 0 | { |
147 | 0 | return GDALDitherRGB2PCTInternal(hRed, hGreen, hBlue, hTarget, hColorTable, |
148 | 0 | 5, nullptr, TRUE, pfnProgress, |
149 | 0 | pProgressArg); |
150 | 0 | } |
151 | | |
152 | | 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 | 0 | { |
159 | 0 | VALIDATE_POINTER1(hRed, "GDALDitherRGB2PCT", CE_Failure); |
160 | 0 | VALIDATE_POINTER1(hGreen, "GDALDitherRGB2PCT", CE_Failure); |
161 | 0 | VALIDATE_POINTER1(hBlue, "GDALDitherRGB2PCT", CE_Failure); |
162 | 0 | VALIDATE_POINTER1(hTarget, "GDALDitherRGB2PCT", CE_Failure); |
163 | 0 | VALIDATE_POINTER1(hColorTable, "GDALDitherRGB2PCT", CE_Failure); |
164 | | |
165 | | /* -------------------------------------------------------------------- */ |
166 | | /* Validate parameters. */ |
167 | | /* -------------------------------------------------------------------- */ |
168 | 0 | const int nXSize = GDALGetRasterBandXSize(hRed); |
169 | 0 | const int nYSize = GDALGetRasterBandYSize(hRed); |
170 | |
|
171 | 0 | if (GDALGetRasterBandXSize(hGreen) != nXSize || |
172 | 0 | GDALGetRasterBandYSize(hGreen) != nYSize || |
173 | 0 | GDALGetRasterBandXSize(hBlue) != nXSize || |
174 | 0 | GDALGetRasterBandYSize(hBlue) != nYSize) |
175 | 0 | { |
176 | 0 | CPLError(CE_Failure, CPLE_IllegalArg, |
177 | 0 | "Green or blue band doesn't match size of red band."); |
178 | |
|
179 | 0 | return CE_Failure; |
180 | 0 | } |
181 | | |
182 | 0 | if (GDALGetRasterBandXSize(hTarget) != nXSize || |
183 | 0 | GDALGetRasterBandYSize(hTarget) != nYSize) |
184 | 0 | { |
185 | 0 | CPLError(CE_Failure, CPLE_IllegalArg, |
186 | 0 | "GDALDitherRGB2PCT(): " |
187 | 0 | "Target band doesn't match size of source bands."); |
188 | |
|
189 | 0 | return CE_Failure; |
190 | 0 | } |
191 | | |
192 | 0 | if (pfnProgress == nullptr) |
193 | 0 | pfnProgress = GDALDummyProgress; |
194 | | |
195 | | /* -------------------------------------------------------------------- */ |
196 | | /* Setup more direct colormap. */ |
197 | | /* -------------------------------------------------------------------- */ |
198 | 0 | int iColor; |
199 | 0 | #ifdef USE_SSE2 |
200 | 0 | int anPCTUnaligned[256 + 4]; // 4 for alignment on 16-byte boundary. |
201 | 0 | int *anPCT = ALIGN_INT_ARRAY_ON_16_BYTE(anPCTUnaligned); |
202 | | #else |
203 | | int anPCT[256 * 4] = {}; |
204 | | #endif |
205 | 0 | const int nColors = GDALGetColorEntryCount(hColorTable); |
206 | |
|
207 | 0 | if (nColors == 0) |
208 | 0 | { |
209 | 0 | CPLError(CE_Failure, CPLE_IllegalArg, |
210 | 0 | "GDALDitherRGB2PCT(): " |
211 | 0 | "Color table must not be empty."); |
212 | |
|
213 | 0 | return CE_Failure; |
214 | 0 | } |
215 | 0 | else if (nColors > 256) |
216 | 0 | { |
217 | 0 | CPLError(CE_Failure, CPLE_IllegalArg, |
218 | 0 | "GDALDitherRGB2PCT(): " |
219 | 0 | "Color table cannot have more than 256 entries."); |
220 | |
|
221 | 0 | return CE_Failure; |
222 | 0 | } |
223 | | |
224 | 0 | iColor = 0; |
225 | 0 | do |
226 | 0 | { |
227 | 0 | GDALColorEntry sEntry; |
228 | |
|
229 | 0 | GDALGetColorEntryAsRGB(hColorTable, iColor, &sEntry); |
230 | 0 | CAST_PCT(anPCT)[4 * iColor + 0] = static_cast<GByte>(sEntry.c1); |
231 | 0 | CAST_PCT(anPCT)[4 * iColor + 1] = static_cast<GByte>(sEntry.c2); |
232 | 0 | CAST_PCT(anPCT)[4 * iColor + 2] = static_cast<GByte>(sEntry.c3); |
233 | 0 | CAST_PCT(anPCT)[4 * iColor + 3] = 0; |
234 | |
|
235 | 0 | iColor++; |
236 | 0 | } while (iColor < nColors); |
237 | |
|
238 | 0 | #ifdef USE_SSE2 |
239 | | // Pad to multiple of 8 colors. |
240 | 0 | const int nColorsMod8 = nColors % 8; |
241 | 0 | if (nColorsMod8) |
242 | 0 | { |
243 | 0 | int iDest = nColors; |
244 | 0 | for (iColor = 0; iColor < 8 - nColorsMod8 && iDest < 256; |
245 | 0 | iColor++, iDest++) |
246 | 0 | { |
247 | 0 | anPCT[iDest] = anPCT[nColors - 1]; |
248 | 0 | } |
249 | 0 | } |
250 | 0 | #endif |
251 | | |
252 | | /* -------------------------------------------------------------------- */ |
253 | | /* Setup various variables. */ |
254 | | /* -------------------------------------------------------------------- */ |
255 | 0 | const int nCLevels = 1 << nBits; |
256 | 0 | const int nCLevelsCube = nCLevels * nCLevels * nCLevels; |
257 | 0 | ColorIndex *psColorIndexMap = nullptr; |
258 | |
|
259 | 0 | GByte *pabyRed = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)); |
260 | 0 | GByte *pabyGreen = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)); |
261 | 0 | GByte *pabyBlue = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)); |
262 | |
|
263 | 0 | GByte *pabyIndex = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize)); |
264 | |
|
265 | 0 | int *panError = |
266 | 0 | static_cast<int *>(VSI_CALLOC_VERBOSE(sizeof(int), (nXSize + 2) * 3)); |
267 | |
|
268 | 0 | if (pabyRed == nullptr || pabyGreen == nullptr || pabyBlue == nullptr || |
269 | 0 | pabyIndex == nullptr || panError == nullptr) |
270 | 0 | { |
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 | 0 | } |
279 | | |
280 | 0 | GByte *pabyColorMap = nullptr; |
281 | 0 | if (pasDynamicColorMap == nullptr) |
282 | 0 | { |
283 | | /* -------------------------------------------------------------------- |
284 | | */ |
285 | | /* Build a 24bit to 8 bit color mapping. */ |
286 | | /* -------------------------------------------------------------------- |
287 | | */ |
288 | |
|
289 | 0 | pabyColorMap = static_cast<GByte *>( |
290 | 0 | VSI_MALLOC_VERBOSE(nCLevelsCube * sizeof(GByte))); |
291 | 0 | if (pabyColorMap == nullptr) |
292 | 0 | { |
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 | 0 | } |
302 | | |
303 | 0 | FindNearestColor(nColors, anPCT, pabyColorMap, nCLevels); |
304 | 0 | } |
305 | 0 | else |
306 | 0 | { |
307 | 0 | pabyColorMap = nullptr; |
308 | 0 | if (nBits == 8 && static_cast<GIntBig>(nXSize) * nYSize <= 65536) |
309 | 0 | { |
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 | 0 | psColorIndexMap = |
314 | 0 | reinterpret_cast<ColorIndex *>(pasDynamicColorMap); |
315 | 0 | memset(psColorIndexMap, 0xFF, sizeof(ColorIndex) * PRIME_FOR_65536); |
316 | 0 | } |
317 | 0 | else |
318 | 0 | { |
319 | 0 | memset(pasDynamicColorMap, 0xFF, 256 * 256 * 256 * sizeof(GInt16)); |
320 | 0 | } |
321 | 0 | } |
322 | | |
323 | | /* ==================================================================== */ |
324 | | /* Loop over all scanlines of data to process. */ |
325 | | /* ==================================================================== */ |
326 | 0 | CPLErr err = CE_None; |
327 | |
|
328 | 0 | for (int iScanline = 0; iScanline < nYSize; iScanline++) |
329 | 0 | { |
330 | | /* -------------------------------------------------------------------- |
331 | | */ |
332 | | /* Report progress */ |
333 | | /* -------------------------------------------------------------------- |
334 | | */ |
335 | 0 | if (!pfnProgress(iScanline / static_cast<double>(nYSize), nullptr, |
336 | 0 | pProgressArg)) |
337 | 0 | { |
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 | 0 | } |
348 | | |
349 | | /* -------------------------------------------------------------------- |
350 | | */ |
351 | | /* Read source data. */ |
352 | | /* -------------------------------------------------------------------- |
353 | | */ |
354 | 0 | CPLErr err1 = GDALRasterIO(hRed, GF_Read, 0, iScanline, nXSize, 1, |
355 | 0 | pabyRed, nXSize, 1, GDT_Byte, 0, 0); |
356 | 0 | if (err1 == CE_None) |
357 | 0 | err1 = GDALRasterIO(hGreen, GF_Read, 0, iScanline, nXSize, 1, |
358 | 0 | pabyGreen, nXSize, 1, GDT_Byte, 0, 0); |
359 | 0 | if (err1 == CE_None) |
360 | 0 | err1 = GDALRasterIO(hBlue, GF_Read, 0, iScanline, nXSize, 1, |
361 | 0 | pabyBlue, nXSize, 1, GDT_Byte, 0, 0); |
362 | 0 | if (err1 != CE_None) |
363 | 0 | { |
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 | 0 | } |
373 | | |
374 | | /* -------------------------------------------------------------------- |
375 | | */ |
376 | | /* Apply the error from the previous line to this one. */ |
377 | | /* -------------------------------------------------------------------- |
378 | | */ |
379 | 0 | if (bDither) |
380 | 0 | { |
381 | 0 | for (int i = 0; i < nXSize; i++) |
382 | 0 | { |
383 | 0 | pabyRed[i] = static_cast<GByte>(std::max( |
384 | 0 | 0, std::min(255, (pabyRed[i] + panError[i * 3 + 0 + 3])))); |
385 | 0 | pabyGreen[i] = static_cast<GByte>(std::max( |
386 | 0 | 0, |
387 | 0 | std::min(255, (pabyGreen[i] + panError[i * 3 + 1 + 3])))); |
388 | 0 | pabyBlue[i] = static_cast<GByte>(std::max( |
389 | 0 | 0, std::min(255, (pabyBlue[i] + panError[i * 3 + 2 + 3])))); |
390 | 0 | } |
391 | |
|
392 | 0 | memset(panError, 0, sizeof(int) * (nXSize + 2) * 3); |
393 | 0 | } |
394 | | |
395 | | /* -------------------------------------------------------------------- |
396 | | */ |
397 | | /* Figure out the nearest color to the RGB value. */ |
398 | | /* -------------------------------------------------------------------- |
399 | | */ |
400 | 0 | int nLastRedError = 0; |
401 | 0 | int nLastGreenError = 0; |
402 | 0 | int nLastBlueError = 0; |
403 | |
|
404 | 0 | for (int i = 0; i < nXSize; i++) |
405 | 0 | { |
406 | 0 | const int nRedValue = |
407 | 0 | std::max(0, std::min(255, pabyRed[i] + nLastRedError)); |
408 | 0 | const int nGreenValue = |
409 | 0 | std::max(0, std::min(255, pabyGreen[i] + nLastGreenError)); |
410 | 0 | const int nBlueValue = |
411 | 0 | std::max(0, std::min(255, pabyBlue[i] + nLastBlueError)); |
412 | |
|
413 | 0 | int iIndex = 0; |
414 | 0 | int nError = 0; |
415 | 0 | int nSixth = 0; |
416 | 0 | if (psColorIndexMap) |
417 | 0 | { |
418 | 0 | const GUInt32 nColorCode = |
419 | 0 | MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue); |
420 | 0 | GUInt32 nIdx = nColorCode % PRIME_FOR_65536; |
421 | 0 | while (true) |
422 | 0 | { |
423 | 0 | if (psColorIndexMap[nIdx].nColorCode == nColorCode) |
424 | 0 | { |
425 | 0 | iIndex = psColorIndexMap[nIdx].nIndex; |
426 | 0 | break; |
427 | 0 | } |
428 | 0 | if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode)) |
429 | 0 | { |
430 | 0 | psColorIndexMap[nIdx].nColorCode = nColorCode; |
431 | 0 | iIndex = FindNearestColor(nColors, anPCT, nRedValue, |
432 | 0 | nGreenValue, nBlueValue); |
433 | 0 | psColorIndexMap[nIdx].nIndex = |
434 | 0 | static_cast<GByte>(iIndex); |
435 | 0 | break; |
436 | 0 | } |
437 | 0 | if (psColorIndexMap[nIdx].nColorCode2 == nColorCode) |
438 | 0 | { |
439 | 0 | iIndex = psColorIndexMap[nIdx].nIndex2; |
440 | 0 | break; |
441 | 0 | } |
442 | 0 | if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode2)) |
443 | 0 | { |
444 | 0 | psColorIndexMap[nIdx].nColorCode2 = nColorCode; |
445 | 0 | iIndex = FindNearestColor(nColors, anPCT, nRedValue, |
446 | 0 | nGreenValue, nBlueValue); |
447 | 0 | psColorIndexMap[nIdx].nIndex2 = |
448 | 0 | static_cast<GByte>(iIndex); |
449 | 0 | break; |
450 | 0 | } |
451 | 0 | if (psColorIndexMap[nIdx].nColorCode3 == nColorCode) |
452 | 0 | { |
453 | 0 | iIndex = psColorIndexMap[nIdx].nIndex3; |
454 | 0 | break; |
455 | 0 | } |
456 | 0 | if (!IsColorCodeSet(psColorIndexMap[nIdx].nColorCode3)) |
457 | 0 | { |
458 | 0 | psColorIndexMap[nIdx].nColorCode3 = nColorCode; |
459 | 0 | iIndex = FindNearestColor(nColors, anPCT, nRedValue, |
460 | 0 | nGreenValue, nBlueValue); |
461 | 0 | psColorIndexMap[nIdx].nIndex3 = |
462 | 0 | static_cast<GByte>(iIndex); |
463 | 0 | break; |
464 | 0 | } |
465 | | |
466 | 0 | do |
467 | 0 | { |
468 | 0 | nIdx += 257; |
469 | 0 | if (nIdx >= PRIME_FOR_65536) |
470 | 0 | nIdx -= PRIME_FOR_65536; |
471 | 0 | } while ( |
472 | 0 | IsColorCodeSet(psColorIndexMap[nIdx].nColorCode) && |
473 | 0 | psColorIndexMap[nIdx].nColorCode != nColorCode && |
474 | 0 | IsColorCodeSet(psColorIndexMap[nIdx].nColorCode2) && |
475 | 0 | psColorIndexMap[nIdx].nColorCode2 != nColorCode && |
476 | 0 | IsColorCodeSet(psColorIndexMap[nIdx].nColorCode3) && |
477 | 0 | psColorIndexMap[nIdx].nColorCode3 != nColorCode); |
478 | 0 | } |
479 | 0 | } |
480 | 0 | else if (pasDynamicColorMap == nullptr) |
481 | 0 | { |
482 | 0 | const int iRed = nRedValue * nCLevels / 256; |
483 | 0 | const int iGreen = nGreenValue * nCLevels / 256; |
484 | 0 | const int iBlue = nBlueValue * nCLevels / 256; |
485 | |
|
486 | 0 | iIndex = pabyColorMap[iRed + iGreen * nCLevels + |
487 | 0 | iBlue * nCLevels * nCLevels]; |
488 | 0 | } |
489 | 0 | else |
490 | 0 | { |
491 | 0 | const GUInt32 nColorCode = |
492 | 0 | MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue); |
493 | 0 | GInt16 *psIndex = &pasDynamicColorMap[nColorCode]; |
494 | 0 | if (*psIndex < 0) |
495 | 0 | { |
496 | 0 | *psIndex = static_cast<GInt16>(FindNearestColor( |
497 | 0 | nColors, anPCT, nRedValue, nGreenValue, nBlueValue)); |
498 | 0 | iIndex = *psIndex; |
499 | 0 | } |
500 | 0 | else |
501 | 0 | { |
502 | 0 | iIndex = *psIndex; |
503 | 0 | } |
504 | 0 | } |
505 | |
|
506 | 0 | pabyIndex[i] = static_cast<GByte>(iIndex); |
507 | 0 | if (!bDither) |
508 | 0 | continue; |
509 | | |
510 | | /* -------------------------------------------------------------------- |
511 | | */ |
512 | | /* Compute Red error, and carry it on to the next error line. |
513 | | */ |
514 | | /* -------------------------------------------------------------------- |
515 | | */ |
516 | 0 | nError = nRedValue - CAST_PCT(anPCT)[4 * iIndex + 0]; |
517 | 0 | nSixth = nError / 6; |
518 | |
|
519 | 0 | panError[i * 3] += nSixth; |
520 | 0 | panError[i * 3 + 6] = nSixth; |
521 | 0 | panError[i * 3 + 3] += nError - 5 * nSixth; |
522 | |
|
523 | 0 | nLastRedError = 2 * nSixth; |
524 | | |
525 | | /* -------------------------------------------------------------------- |
526 | | */ |
527 | | /* Compute Green error, and carry it on to the next error line. |
528 | | */ |
529 | | /* -------------------------------------------------------------------- |
530 | | */ |
531 | 0 | nError = nGreenValue - CAST_PCT(anPCT)[4 * iIndex + 1]; |
532 | 0 | nSixth = nError / 6; |
533 | |
|
534 | 0 | panError[i * 3 + 1] += nSixth; |
535 | 0 | panError[i * 3 + 6 + 1] = nSixth; |
536 | 0 | panError[i * 3 + 3 + 1] += nError - 5 * nSixth; |
537 | |
|
538 | 0 | nLastGreenError = 2 * nSixth; |
539 | | |
540 | | /* -------------------------------------------------------------------- |
541 | | */ |
542 | | /* Compute Blue error, and carry it on to the next error line. |
543 | | */ |
544 | | /* -------------------------------------------------------------------- |
545 | | */ |
546 | 0 | nError = nBlueValue - CAST_PCT(anPCT)[4 * iIndex + 2]; |
547 | 0 | nSixth = nError / 6; |
548 | |
|
549 | 0 | panError[i * 3 + 2] += nSixth; |
550 | 0 | panError[i * 3 + 6 + 2] = nSixth; |
551 | 0 | panError[i * 3 + 3 + 2] += nError - 5 * nSixth; |
552 | |
|
553 | 0 | nLastBlueError = 2 * nSixth; |
554 | 0 | } |
555 | | |
556 | | /* -------------------------------------------------------------------- |
557 | | */ |
558 | | /* Write results. */ |
559 | | /* -------------------------------------------------------------------- |
560 | | */ |
561 | 0 | err = GDALRasterIO(hTarget, GF_Write, 0, iScanline, nXSize, 1, |
562 | 0 | pabyIndex, nXSize, 1, GDT_Byte, 0, 0); |
563 | 0 | if (err != CE_None) |
564 | 0 | break; |
565 | 0 | } |
566 | | |
567 | 0 | pfnProgress(1.0, nullptr, pProgressArg); |
568 | | |
569 | | /* -------------------------------------------------------------------- */ |
570 | | /* Cleanup */ |
571 | | /* -------------------------------------------------------------------- */ |
572 | 0 | CPLFree(pabyRed); |
573 | 0 | CPLFree(pabyGreen); |
574 | 0 | CPLFree(pabyBlue); |
575 | 0 | CPLFree(pabyIndex); |
576 | 0 | CPLFree(panError); |
577 | 0 | CPLFree(pabyColorMap); |
578 | |
|
579 | 0 | return err; |
580 | 0 | } |
581 | | |
582 | | static int FindNearestColor(int nColors, int *panPCT, int nRedValue, |
583 | | int nGreenValue, int nBlueValue) |
584 | | |
585 | 0 | { |
586 | 0 | #ifdef USE_SSE2 |
587 | 0 | int nBestDist = 768; |
588 | 0 | int nBestIndex = 0; |
589 | |
|
590 | 0 | int anDistanceUnaligned[16 + 4] = |
591 | 0 | {}; // 4 for alignment on 16-byte boundary. |
592 | 0 | int *anDistance = ALIGN_INT_ARRAY_ON_16_BYTE(anDistanceUnaligned); |
593 | |
|
594 | 0 | const __m128i ff = _mm_set1_epi32(0xFFFFFFFF); |
595 | 0 | const __m128i mask_low = _mm_srli_epi64(ff, 32); |
596 | 0 | const __m128i mask_high = _mm_slli_epi64(ff, 32); |
597 | |
|
598 | 0 | const unsigned int nColorVal = |
599 | 0 | MAKE_COLOR_CODE(nRedValue, nGreenValue, nBlueValue); |
600 | 0 | const __m128i thisColor = _mm_set1_epi32(nColorVal); |
601 | 0 | const __m128i thisColor_low = _mm_srli_epi64(thisColor, 32); |
602 | 0 | const __m128i thisColor_high = _mm_slli_epi64(thisColor, 32); |
603 | |
|
604 | 0 | for (int iColor = 0; iColor < nColors; iColor += 8) |
605 | 0 | { |
606 | 0 | const __m128i pctColor = |
607 | 0 | _mm_load_si128(reinterpret_cast<__m128i *>(&panPCT[iColor])); |
608 | 0 | const __m128i pctColor2 = |
609 | 0 | _mm_load_si128(reinterpret_cast<__m128i *>(&panPCT[iColor + 4])); |
610 | |
|
611 | 0 | _mm_store_si128( |
612 | 0 | reinterpret_cast<__m128i *>(anDistance), |
613 | 0 | _mm_sad_epu8(_mm_and_si128(pctColor, mask_low), thisColor_low)); |
614 | 0 | _mm_store_si128( |
615 | 0 | reinterpret_cast<__m128i *>(anDistance + 4), |
616 | 0 | _mm_sad_epu8(_mm_and_si128(pctColor, mask_high), thisColor_high)); |
617 | 0 | _mm_store_si128( |
618 | 0 | reinterpret_cast<__m128i *>(anDistance + 8), |
619 | 0 | _mm_sad_epu8(_mm_and_si128(pctColor2, mask_low), thisColor_low)); |
620 | 0 | _mm_store_si128( |
621 | 0 | reinterpret_cast<__m128i *>(anDistance + 12), |
622 | 0 | _mm_sad_epu8(_mm_and_si128(pctColor2, mask_high), thisColor_high)); |
623 | |
|
624 | 0 | if (anDistance[0] < nBestDist) |
625 | 0 | { |
626 | 0 | nBestIndex = iColor; |
627 | 0 | nBestDist = anDistance[0]; |
628 | 0 | } |
629 | 0 | if (anDistance[4] < nBestDist) |
630 | 0 | { |
631 | 0 | nBestIndex = iColor + 1; |
632 | 0 | nBestDist = anDistance[4]; |
633 | 0 | } |
634 | 0 | if (anDistance[2] < nBestDist) |
635 | 0 | { |
636 | 0 | nBestIndex = iColor + 2; |
637 | 0 | nBestDist = anDistance[2]; |
638 | 0 | } |
639 | 0 | if (anDistance[6] < nBestDist) |
640 | 0 | { |
641 | 0 | nBestIndex = iColor + 3; |
642 | 0 | nBestDist = anDistance[6]; |
643 | 0 | } |
644 | 0 | if (anDistance[8 + 0] < nBestDist) |
645 | 0 | { |
646 | 0 | nBestIndex = iColor + 4; |
647 | 0 | nBestDist = anDistance[8 + 0]; |
648 | 0 | } |
649 | 0 | if (anDistance[8 + 4] < nBestDist) |
650 | 0 | { |
651 | 0 | nBestIndex = iColor + 4 + 1; |
652 | 0 | nBestDist = anDistance[8 + 4]; |
653 | 0 | } |
654 | 0 | if (anDistance[8 + 2] < nBestDist) |
655 | 0 | { |
656 | 0 | nBestIndex = iColor + 4 + 2; |
657 | 0 | nBestDist = anDistance[8 + 2]; |
658 | 0 | } |
659 | 0 | if (anDistance[8 + 6] < nBestDist) |
660 | 0 | { |
661 | 0 | nBestIndex = iColor + 4 + 3; |
662 | 0 | nBestDist = anDistance[8 + 6]; |
663 | 0 | } |
664 | 0 | } |
665 | 0 | 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 | 0 | } |
685 | | |
686 | | /************************************************************************/ |
687 | | /* FindNearestColor() */ |
688 | | /* */ |
689 | | /* Finear near PCT color for any RGB color. */ |
690 | | /************************************************************************/ |
691 | | |
692 | | static void FindNearestColor(int nColors, int *panPCT, GByte *pabyColorMap, |
693 | | int nCLevels) |
694 | | |
695 | 0 | { |
696 | | /* -------------------------------------------------------------------- */ |
697 | | /* Loop over all the cells in the high density cube. */ |
698 | | /* -------------------------------------------------------------------- */ |
699 | 0 | for (int iBlue = 0; iBlue < nCLevels; iBlue++) |
700 | 0 | { |
701 | 0 | for (int iGreen = 0; iGreen < nCLevels; iGreen++) |
702 | 0 | { |
703 | 0 | for (int iRed = 0; iRed < nCLevels; iRed++) |
704 | 0 | { |
705 | 0 | const int nRedValue = (iRed * 255) / (nCLevels - 1); |
706 | 0 | const int nGreenValue = (iGreen * 255) / (nCLevels - 1); |
707 | 0 | const int nBlueValue = (iBlue * 255) / (nCLevels - 1); |
708 | |
|
709 | 0 | const int nBestIndex = FindNearestColor( |
710 | 0 | nColors, panPCT, nRedValue, nGreenValue, nBlueValue); |
711 | 0 | pabyColorMap[iRed + iGreen * nCLevels + |
712 | 0 | iBlue * nCLevels * nCLevels] = |
713 | 0 | static_cast<GByte>(nBestIndex); |
714 | 0 | } |
715 | 0 | } |
716 | 0 | } |
717 | 0 | } |