Coverage Report

Created: 2025-06-09 07:29

/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
}