Coverage Report

Created: 2025-06-09 07:38

/src/gdal/alg/llrasterize.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  GDAL
4
 * Purpose:  Vector polygon rasterization code.
5
 * Author:   Frank Warmerdam, warmerdam@pobox.com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
9
 * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 ****************************************************************************/
13
14
#include "cpl_port.h"
15
#include "gdal_alg_priv.h"
16
17
#include <cmath>
18
#include <cstdlib>
19
#include <cstring>
20
21
#include <algorithm>
22
#include <set>
23
#include <utility>
24
#include <vector>
25
26
#include "gdal_alg.h"
27
28
/************************************************************************/
29
/*                       dllImageFilledPolygon()                        */
30
/*                                                                      */
31
/*      Perform scanline conversion of the passed multi-ring            */
32
/*      polygon.  Note the polygon does not need to be explicitly       */
33
/*      closed.  The scanline function will be called with              */
34
/*      horizontal scanline chunks which may not be entirely            */
35
/*      contained within the valid raster area (in the X                */
36
/*      direction).                                                     */
37
/*                                                                      */
38
/*      NEW: Nodes' coordinate are kept as double  in order             */
39
/*      to compute accurately the intersections with the lines          */
40
/*                                                                      */
41
/*      A pixel is considered inside a polygon if its center            */
42
/*      falls inside the polygon. This is robust unless                 */
43
/*      the nodes are placed in the center of the pixels in which       */
44
/*      case, due to numerical inaccuracies, it's hard to predict       */
45
/*      if the pixel will be considered inside or outside the shape.    */
46
/************************************************************************/
47
48
/*
49
 * NOTE: This code was originally adapted from the gdImageFilledPolygon()
50
 * function in libgd.
51
 *
52
 * http://www.boutell.com/gd/
53
 *
54
 * It was later adapted for direct inclusion in GDAL and relicensed under
55
 * the GDAL MIT license (pulled from the OpenEV distribution).
56
 */
57
58
void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
59
                               int nPartCount, const int *panPartSize,
60
                               const double *padfX, const double *padfY,
61
                               const double *dfVariant,
62
                               llScanlineFunc pfnScanlineFunc,
63
                               GDALRasterizeInfo *pCBData,
64
                               bool bAvoidBurningSamePoints)
65
0
{
66
0
    if (!nPartCount)
67
0
    {
68
0
        return;
69
0
    }
70
71
0
    int n = 0;
72
0
    for (int part = 0; part < nPartCount; part++)
73
0
        n += panPartSize[part];
74
75
0
    std::vector<int> polyInts(n);
76
0
    std::vector<int> polyInts2;
77
0
    if (bAvoidBurningSamePoints)
78
0
        polyInts2.resize(n);
79
80
0
    double dminy = padfY[0];
81
0
    double dmaxy = padfY[0];
82
0
    for (int i = 1; i < n; i++)
83
0
    {
84
0
        if (padfY[i] < dminy)
85
0
        {
86
0
            dminy = padfY[i];
87
0
        }
88
0
        if (padfY[i] > dmaxy)
89
0
        {
90
0
            dmaxy = padfY[i];
91
0
        }
92
0
    }
93
0
    int miny = static_cast<int>(dminy);
94
0
    int maxy = static_cast<int>(dmaxy);
95
96
0
    if (miny < 0)
97
0
        miny = 0;
98
0
    if (maxy >= nRasterYSize)
99
0
        maxy = nRasterYSize - 1;
100
101
0
    int minx = 0;
102
0
    const int maxx = nRasterXSize - 1;
103
104
    // Fix in 1.3: count a vertex only once.
105
0
    for (int y = miny; y <= maxy; y++)
106
0
    {
107
0
        int partoffset = 0;
108
109
0
        const double dy = y + 0.5;  // Center height of line.
110
111
0
        int part = 0;
112
0
        int ints = 0;
113
0
        int ints2 = 0;
114
115
0
        for (int i = 0; i < n; i++)
116
0
        {
117
0
            if (i == partoffset + panPartSize[part])
118
0
            {
119
0
                partoffset += panPartSize[part];
120
0
                part++;
121
0
            }
122
123
0
            int ind1 = 0;
124
0
            int ind2 = 0;
125
0
            if (i == partoffset)
126
0
            {
127
0
                ind1 = partoffset + panPartSize[part] - 1;
128
0
                ind2 = partoffset;
129
0
            }
130
0
            else
131
0
            {
132
0
                ind1 = i - 1;
133
0
                ind2 = i;
134
0
            }
135
136
0
            double dy1 = padfY[ind1];
137
0
            double dy2 = padfY[ind2];
138
139
0
            if ((dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy))
140
0
                continue;
141
142
0
            double dx1 = 0.0;
143
0
            double dx2 = 0.0;
144
0
            if (dy1 < dy2)
145
0
            {
146
0
                dx1 = padfX[ind1];
147
0
                dx2 = padfX[ind2];
148
0
            }
149
0
            else if (dy1 > dy2)
150
0
            {
151
0
                std::swap(dy1, dy2);
152
0
                dx2 = padfX[ind1];
153
0
                dx1 = padfX[ind2];
154
0
            }
155
0
            else  // if( fabs(dy1-dy2) < 1.e-6 )
156
0
            {
157
158
                // AE: DO NOT skip bottom horizontal segments
159
                // -Fill them separately-
160
0
                if (padfX[ind1] > padfX[ind2])
161
0
                {
162
0
                    const int horizontal_x1 =
163
0
                        static_cast<int>(floor(padfX[ind2] + 0.5));
164
0
                    const int horizontal_x2 =
165
0
                        static_cast<int>(floor(padfX[ind1] + 0.5));
166
167
0
                    if ((horizontal_x1 > maxx) || (horizontal_x2 <= minx))
168
0
                        continue;
169
170
                    // Fill the horizontal segment (separately from the rest).
171
0
                    if (bAvoidBurningSamePoints)
172
0
                    {
173
0
                        polyInts2[ints2++] = horizontal_x1;
174
0
                        polyInts2[ints2++] = horizontal_x2;
175
0
                    }
176
0
                    else
177
0
                    {
178
0
                        pfnScanlineFunc(
179
0
                            pCBData, y, horizontal_x1, horizontal_x2 - 1,
180
0
                            dfVariant == nullptr ? 0 : dfVariant[0]);
181
0
                    }
182
0
                }
183
                // else: Skip top horizontal segments.
184
                // They are already filled in the regular loop.
185
0
                continue;
186
0
            }
187
188
0
            if (dy < dy2 && dy >= dy1)
189
0
            {
190
0
                const double intersect =
191
0
                    (dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1;
192
193
0
                polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
194
0
            }
195
0
        }
196
197
0
        std::sort(polyInts.begin(), polyInts.begin() + ints);
198
0
        std::sort(polyInts2.begin(), polyInts2.begin() + ints2);
199
200
0
        for (int i = 0; i + 1 < ints; i += 2)
201
0
        {
202
0
            if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
203
0
            {
204
0
                pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
205
0
                                dfVariant == nullptr ? 0 : dfVariant[0]);
206
0
            }
207
0
        }
208
209
0
        for (int i2 = 0, i = 0; i2 + 1 < ints2; i2 += 2)
210
0
        {
211
0
            if (polyInts2[i2] <= maxx && polyInts2[i2 + 1] > minx)
212
0
            {
213
                // "synchronize" polyInts[i] with polyInts2[i2]
214
0
                while (i + 1 < ints && polyInts[i] < polyInts2[i2])
215
0
                    i += 2;
216
                // Only burn if we don't have a common segment between
217
                // polyInts[] and polyInts2[]
218
0
                if (i + 1 >= ints || polyInts[i] != polyInts2[i2])
219
0
                {
220
0
                    pfnScanlineFunc(pCBData, y, polyInts2[i2],
221
0
                                    polyInts2[i2 + 1] - 1,
222
0
                                    dfVariant == nullptr ? 0 : dfVariant[0]);
223
0
                }
224
0
            }
225
0
        }
226
0
    }
227
0
}
228
229
/************************************************************************/
230
/*                         GDALdllImagePoint()                          */
231
/************************************************************************/
232
233
void GDALdllImagePoint(int nRasterXSize, int nRasterYSize, int nPartCount,
234
                       const int * /*panPartSize*/, const double *padfX,
235
                       const double *padfY, const double *padfVariant,
236
                       llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
237
0
{
238
0
    for (int i = 0; i < nPartCount; i++)
239
0
    {
240
0
        const int nX = static_cast<int>(floor(padfX[i]));
241
0
        const int nY = static_cast<int>(floor(padfY[i]));
242
0
        double dfVariant = 0.0;
243
0
        if (padfVariant != nullptr)
244
0
            dfVariant = padfVariant[i];
245
246
0
        if (0 <= nX && nX < nRasterXSize && 0 <= nY && nY < nRasterYSize)
247
0
            pfnPointFunc(pCBData, nY, nX, dfVariant);
248
0
    }
249
0
}
250
251
/************************************************************************/
252
/*                         GDALdllImageLine()                           */
253
/************************************************************************/
254
255
void GDALdllImageLine(int nRasterXSize, int nRasterYSize, int nPartCount,
256
                      const int *panPartSize, const double *padfX,
257
                      const double *padfY, const double *padfVariant,
258
                      llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
259
0
{
260
0
    if (!nPartCount)
261
0
        return;
262
263
0
    for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
264
0
    {
265
0
        for (int j = 1; j < panPartSize[i]; j++)
266
0
        {
267
0
            int iX = static_cast<int>(floor(padfX[n + j - 1]));
268
0
            int iY = static_cast<int>(floor(padfY[n + j - 1]));
269
270
0
            const int iX1 = static_cast<int>(floor(padfX[n + j]));
271
0
            const int iY1 = static_cast<int>(floor(padfY[n + j]));
272
273
0
            double dfVariant = 0.0;
274
0
            double dfVariant1 = 0.0;
275
0
            if (padfVariant != nullptr &&
276
0
                pCBData->eBurnValueSource != GBV_UserBurnValue)
277
0
            {
278
0
                dfVariant = padfVariant[n + j - 1];
279
0
                dfVariant1 = padfVariant[n + j];
280
0
            }
281
282
0
            int nDeltaX = std::abs(iX1 - iX);
283
0
            int nDeltaY = std::abs(iY1 - iY);
284
285
            // Step direction depends on line direction.
286
0
            const int nXStep = (iX > iX1) ? -1 : 1;
287
0
            const int nYStep = (iY > iY1) ? -1 : 1;
288
289
            // Determine the line slope.
290
0
            if (nDeltaX >= nDeltaY)
291
0
            {
292
0
                const int nXError = nDeltaY << 1;
293
0
                const int nYError = nXError - (nDeltaX << 1);
294
0
                int nError = nXError - nDeltaX;
295
                // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
296
                // but if it is == 0, dfDeltaVariant is not really used, so any
297
                // value is okay.
298
0
                const double dfDeltaVariant =
299
0
                    nDeltaX == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaX;
300
301
                // Do not burn the end point, unless we are in the last
302
                // segment. This is to avoid burning twice intermediate points,
303
                // which causes artifacts in Add mode
304
0
                if (j != panPartSize[i] - 1)
305
0
                {
306
0
                    nDeltaX--;
307
0
                }
308
309
0
                while (nDeltaX-- >= 0)
310
0
                {
311
0
                    if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
312
0
                        iY < nRasterYSize)
313
0
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
314
315
0
                    dfVariant += dfDeltaVariant;
316
0
                    iX += nXStep;
317
0
                    if (nError > 0)
318
0
                    {
319
0
                        iY += nYStep;
320
0
                        nError += nYError;
321
0
                    }
322
0
                    else
323
0
                    {
324
0
                        nError += nXError;
325
0
                    }
326
0
                }
327
0
            }
328
0
            else
329
0
            {
330
0
                const int nXError = nDeltaX << 1;
331
0
                const int nYError = nXError - (nDeltaY << 1);
332
0
                int nError = nXError - nDeltaY;
333
                // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
334
                // but if it is == 0, dfDeltaVariant is not really used, so any
335
                // value is okay.
336
0
                double dfDeltaVariant =
337
0
                    nDeltaY == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaY;
338
339
                // Do not burn the end point, unless we are in the last
340
                // segment. This is to avoid burning twice intermediate points,
341
                // which causes artifacts in Add mode
342
0
                if (j != panPartSize[i] - 1)
343
0
                {
344
0
                    nDeltaY--;
345
0
                }
346
347
0
                while (nDeltaY-- >= 0)
348
0
                {
349
0
                    if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
350
0
                        iY < nRasterYSize)
351
0
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
352
353
0
                    dfVariant += dfDeltaVariant;
354
0
                    iY += nYStep;
355
0
                    if (nError > 0)
356
0
                    {
357
0
                        iX += nXStep;
358
0
                        nError += nYError;
359
0
                    }
360
0
                    else
361
0
                    {
362
0
                        nError += nXError;
363
0
                    }
364
0
                }
365
0
            }
366
0
        }
367
0
    }
368
0
}
369
370
/************************************************************************/
371
/*                     GDALdllImageLineAllTouched()                     */
372
/*                                                                      */
373
/*      This alternate line drawing algorithm attempts to ensure        */
374
/*      that every pixel touched at all by the line will get set.       */
375
/*      @param padfVariant should contain the values that are to be     */
376
/*      added to the burn value.  The values along the line between the */
377
/*      points will be linearly interpolated. These values are used     */
378
/*      only if pCBData->eBurnValueSource is set to something other     */
379
/*      than GBV_UserBurnValue. If NULL is passed, a monotonous line    */
380
/*      will be drawn with the burn value.                              */
381
/************************************************************************/
382
383
void GDALdllImageLineAllTouched(
384
    int nRasterXSize, int nRasterYSize, int nPartCount, const int *panPartSize,
385
    const double *padfX, const double *padfY, const double *padfVariant,
386
    llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData,
387
    bool bAvoidBurningSamePoints, bool bIntersectOnly)
388
389
0
{
390
    // This is an epsilon to detect geometries that are aligned with pixel
391
    // coordinates. Hard to find the right value. We put it to that value
392
    // to satisfy the scenarios of https://github.com/OSGeo/gdal/issues/7523
393
    // and https://github.com/OSGeo/gdal/issues/6414
394
0
    constexpr double EPSILON_INTERSECT_ONLY = 1e-4;
395
396
0
    if (!nPartCount)
397
0
        return;
398
399
0
    for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
400
0
    {
401
0
        std::set<std::pair<int, int>> lastBurntPoints;
402
0
        std::set<std::pair<int, int>> newBurntPoints;
403
404
0
        for (int j = 1; j < panPartSize[i]; j++)
405
0
        {
406
0
            lastBurntPoints = std::move(newBurntPoints);
407
0
            newBurntPoints.clear();
408
409
0
            double dfX = padfX[n + j - 1];
410
0
            double dfY = padfY[n + j - 1];
411
412
0
            double dfXEnd = padfX[n + j];
413
0
            double dfYEnd = padfY[n + j];
414
415
0
            double dfVariant = 0.0;
416
0
            double dfVariantEnd = 0.0;
417
0
            if (padfVariant != nullptr &&
418
0
                static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
419
0
                    GBV_UserBurnValue)
420
0
            {
421
0
                dfVariant = padfVariant[n + j - 1];
422
0
                dfVariantEnd = padfVariant[n + j];
423
0
            }
424
425
            // Skip segments that are off the target region.
426
0
            if ((dfY < 0.0 && dfYEnd < 0.0) ||
427
0
                (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
428
0
                (dfX < 0.0 && dfXEnd < 0.0) ||
429
0
                (dfX > nRasterXSize && dfXEnd > nRasterXSize))
430
0
                continue;
431
432
            // Swap if needed so we can proceed from left2right (X increasing)
433
0
            if (dfX > dfXEnd)
434
0
            {
435
0
                std::swap(dfX, dfXEnd);
436
0
                std::swap(dfY, dfYEnd);
437
0
                std::swap(dfVariant, dfVariantEnd);
438
0
            }
439
440
            // Special case for vertical lines.
441
442
0
            if (fabs(dfX - dfXEnd) < .01)
443
0
            {
444
0
                if (bIntersectOnly)
445
0
                {
446
0
                    if (std::abs(dfX - std::round(dfX)) <
447
0
                            EPSILON_INTERSECT_ONLY &&
448
0
                        std::abs(dfXEnd - std::round(dfXEnd)) <
449
0
                            EPSILON_INTERSECT_ONLY)
450
0
                        continue;
451
0
                }
452
453
0
                if (dfYEnd < dfY)
454
0
                {
455
0
                    std::swap(dfY, dfYEnd);
456
0
                    std::swap(dfVariant, dfVariantEnd);
457
0
                }
458
459
0
                const int iX = static_cast<int>(floor(dfXEnd));
460
0
                int iY = static_cast<int>(floor(dfY));
461
0
                int iYEnd =
462
0
                    static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));
463
464
0
                if (iX < 0 || iX >= nRasterXSize)
465
0
                    continue;
466
467
0
                double dfDeltaVariant = 0.0;
468
0
                if (dfYEnd - dfY > 0.0)
469
0
                    dfDeltaVariant = (dfVariantEnd - dfVariant) /
470
0
                                     (dfYEnd - dfY);  // Per unit change in iY.
471
472
                // Clip to the borders of the target region.
473
0
                if (iY < 0)
474
0
                    iY = 0;
475
0
                if (iYEnd >= nRasterYSize)
476
0
                    iYEnd = nRasterYSize - 1;
477
0
                dfVariant += dfDeltaVariant * (iY - dfY);
478
479
0
                if (padfVariant == nullptr)
480
0
                {
481
0
                    for (; iY <= iYEnd; iY++)
482
0
                    {
483
0
                        if (bAvoidBurningSamePoints)
484
0
                        {
485
0
                            auto yx = std::pair<int, int>(iY, iX);
486
0
                            if (lastBurntPoints.find(yx) !=
487
0
                                lastBurntPoints.end())
488
0
                            {
489
0
                                continue;
490
0
                            }
491
0
                            newBurntPoints.insert(yx);
492
0
                        }
493
0
                        pfnPointFunc(pCBData, iY, iX, 0.0);
494
0
                    }
495
0
                }
496
0
                else
497
0
                {
498
0
                    for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
499
0
                    {
500
0
                        if (bAvoidBurningSamePoints)
501
0
                        {
502
0
                            auto yx = std::pair<int, int>(iY, iX);
503
0
                            if (lastBurntPoints.find(yx) !=
504
0
                                lastBurntPoints.end())
505
0
                            {
506
0
                                continue;
507
0
                            }
508
0
                            newBurntPoints.insert(yx);
509
0
                        }
510
0
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
511
0
                    }
512
0
                }
513
514
0
                continue;  // Next segment.
515
0
            }
516
517
0
            const double dfDeltaVariant =
518
0
                (dfVariantEnd - dfVariant) /
519
0
                (dfXEnd - dfX);  // Per unit change in iX.
520
521
            // Special case for horizontal lines.
522
0
            if (fabs(dfY - dfYEnd) < .01)
523
0
            {
524
0
                if (bIntersectOnly)
525
0
                {
526
0
                    if (std::abs(dfY - std::round(dfY)) <
527
0
                            EPSILON_INTERSECT_ONLY &&
528
0
                        std::abs(dfYEnd - std::round(dfYEnd)) <
529
0
                            EPSILON_INTERSECT_ONLY)
530
0
                        continue;
531
0
                }
532
533
0
                if (dfXEnd < dfX)
534
0
                {
535
0
                    std::swap(dfX, dfXEnd);
536
0
                    std::swap(dfVariant, dfVariantEnd);
537
0
                }
538
539
0
                int iX = static_cast<int>(floor(dfX));
540
0
                const int iY = static_cast<int>(floor(dfY));
541
0
                int iXEnd =
542
0
                    static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));
543
544
0
                if (iY < 0 || iY >= nRasterYSize)
545
0
                    continue;
546
547
                // Clip to the borders of the target region.
548
0
                if (iX < 0)
549
0
                    iX = 0;
550
0
                if (iXEnd >= nRasterXSize)
551
0
                    iXEnd = nRasterXSize - 1;
552
0
                dfVariant += dfDeltaVariant * (iX - dfX);
553
554
0
                if (padfVariant == nullptr)
555
0
                {
556
0
                    for (; iX <= iXEnd; iX++)
557
0
                    {
558
0
                        if (bAvoidBurningSamePoints)
559
0
                        {
560
0
                            auto yx = std::pair<int, int>(iY, iX);
561
0
                            if (lastBurntPoints.find(yx) !=
562
0
                                lastBurntPoints.end())
563
0
                            {
564
0
                                continue;
565
0
                            }
566
0
                            newBurntPoints.insert(yx);
567
0
                        }
568
0
                        pfnPointFunc(pCBData, iY, iX, 0.0);
569
0
                    }
570
0
                }
571
0
                else
572
0
                {
573
0
                    for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
574
0
                    {
575
0
                        if (bAvoidBurningSamePoints)
576
0
                        {
577
0
                            auto yx = std::pair<int, int>(iY, iX);
578
0
                            if (lastBurntPoints.find(yx) !=
579
0
                                lastBurntPoints.end())
580
0
                            {
581
0
                                continue;
582
0
                            }
583
0
                            newBurntPoints.insert(yx);
584
0
                        }
585
0
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
586
0
                    }
587
0
                }
588
589
0
                continue;  // Next segment.
590
0
            }
591
592
            /* --------------------------------------------------------------------
593
             */
594
            /*      General case - left to right sloped. */
595
            /* --------------------------------------------------------------------
596
             */
597
0
            const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
598
599
            // Clip segment in X.
600
0
            if (dfXEnd > nRasterXSize)
601
0
            {
602
0
                dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
603
0
                dfXEnd = nRasterXSize;
604
0
            }
605
0
            if (dfX < 0.0)
606
0
            {
607
0
                dfY += (0.0 - dfX) * dfSlope;
608
0
                dfVariant += dfDeltaVariant * (0.0 - dfX);
609
0
                dfX = 0.0;
610
0
            }
611
612
            // Clip segment in Y.
613
0
            if (dfYEnd > dfY)
614
0
            {
615
0
                if (dfY < 0.0)
616
0
                {
617
0
                    const double dfDiffX = (0.0 - dfY) / dfSlope;
618
0
                    dfX += dfDiffX;
619
0
                    dfVariant += dfDeltaVariant * dfDiffX;
620
0
                    dfY = 0.0;
621
0
                }
622
0
                if (dfYEnd >= nRasterYSize)
623
0
                {
624
0
                    dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
625
0
                    if (dfXEnd > nRasterXSize)
626
0
                        dfXEnd = nRasterXSize;
627
                    // dfYEnd is no longer used afterwards, but for
628
                    // consistency it should be:
629
                    // dfYEnd = nRasterXSize;
630
0
                }
631
0
            }
632
0
            else
633
0
            {
634
0
                if (dfY >= nRasterYSize)
635
0
                {
636
0
                    const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
637
0
                    dfX += dfDiffX;
638
0
                    dfVariant += dfDeltaVariant * dfDiffX;
639
0
                    dfY = nRasterYSize;
640
0
                }
641
0
                if (dfYEnd < 0.0)
642
0
                {
643
0
                    dfXEnd -= (dfYEnd - 0) / dfSlope;
644
                    // dfYEnd is no longer used afterwards, but for
645
                    // consistency it should be:
646
                    // dfYEnd = 0.0;
647
0
                }
648
0
            }
649
650
            // Step from pixel to pixel.
651
0
            while (dfX >= 0.0 && dfX < dfXEnd)
652
0
            {
653
0
                const int iX = static_cast<int>(floor(dfX));
654
0
                const int iY = static_cast<int>(floor(dfY));
655
656
                // Burn in the current point.
657
                // We should be able to drop the Y check because we clipped
658
                // in Y, but there may be some error with all the small steps.
659
0
                if (iY >= 0 && iY < nRasterYSize)
660
0
                {
661
0
                    if (bAvoidBurningSamePoints)
662
0
                    {
663
0
                        auto yx = std::pair<int, int>(iY, iX);
664
0
                        if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
665
0
                            newBurntPoints.find(yx) == newBurntPoints.end())
666
0
                        {
667
0
                            newBurntPoints.insert(yx);
668
0
                            pfnPointFunc(pCBData, iY, iX, dfVariant);
669
0
                        }
670
0
                    }
671
0
                    else
672
0
                    {
673
0
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
674
0
                    }
675
0
                }
676
677
0
                double dfStepX = floor(dfX + 1.0) - dfX;
678
0
                double dfStepY = dfStepX * dfSlope;
679
680
                // Step to right pixel without changing scanline?
681
0
                if (static_cast<int>(floor(dfY + dfStepY)) == iY)
682
0
                {
683
0
                    dfX += dfStepX;
684
0
                    dfY += dfStepY;
685
0
                    dfVariant += dfDeltaVariant * dfStepX;
686
0
                }
687
0
                else if (dfSlope < 0)
688
0
                {
689
0
                    dfStepY = iY - dfY;
690
0
                    if (dfStepY > -0.000000001)
691
0
                        dfStepY = -0.000000001;
692
693
0
                    dfStepX = dfStepY / dfSlope;
694
0
                    dfX += dfStepX;
695
0
                    dfY += dfStepY;
696
0
                    dfVariant += dfDeltaVariant * dfStepX;
697
0
                }
698
0
                else
699
0
                {
700
0
                    dfStepY = (iY + 1) - dfY;
701
0
                    if (dfStepY < 0.000000001)
702
0
                        dfStepY = 0.000000001;
703
704
0
                    dfStepX = dfStepY / dfSlope;
705
0
                    dfX += dfStepX;
706
0
                    dfY += dfStepY;
707
0
                    dfVariant += dfDeltaVariant * dfStepX;
708
0
                }
709
0
            }  // Next step along segment.
710
0
        }      // Next segment.
711
0
    }          // Next part.
712
0
}