Line data Source code
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 326 : 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 : {
66 326 : if (!nPartCount)
67 : {
68 0 : return;
69 : }
70 :
71 326 : int n = 0;
72 663 : for (int part = 0; part < nPartCount; part++)
73 337 : n += panPartSize[part];
74 :
75 652 : std::vector<int> polyInts(n);
76 652 : std::vector<int> polyInts2;
77 326 : if (bAvoidBurningSamePoints)
78 15 : polyInts2.resize(n);
79 :
80 326 : double dminy = padfY[0];
81 326 : double dmaxy = padfY[0];
82 9412 : for (int i = 1; i < n; i++)
83 : {
84 9086 : if (padfY[i] < dminy)
85 : {
86 1276 : dminy = padfY[i];
87 : }
88 7810 : else if (padfY[i] > dmaxy)
89 : {
90 2063 : dmaxy = padfY[i];
91 : }
92 : }
93 326 : const int miny = static_cast<int>(std::max(0.0, dminy));
94 : const int maxy =
95 326 : static_cast<int>(std::min<double>(dmaxy, nRasterYSize - 1));
96 :
97 326 : constexpr int minx = 0;
98 326 : const int maxx = nRasterXSize - 1;
99 :
100 : // Fix in 1.3: count a vertex only once.
101 12256 : for (int y = miny; y <= maxy; y++)
102 : {
103 11930 : int partoffset = 0;
104 :
105 11930 : const double dy = y + 0.5; // Center height of line.
106 :
107 11930 : int part = 0;
108 11930 : int ints = 0;
109 11930 : int ints2 = 0;
110 :
111 7037230 : for (int i = 0; i < n; i++)
112 : {
113 7025300 : if (i == partoffset + panPartSize[part])
114 : {
115 150 : partoffset += panPartSize[part];
116 150 : part++;
117 : }
118 :
119 7025300 : int ind1 = 0;
120 7025300 : int ind2 = 0;
121 7025300 : if (i == partoffset)
122 : {
123 12080 : ind1 = partoffset + panPartSize[part] - 1;
124 12080 : ind2 = partoffset;
125 : }
126 : else
127 : {
128 7013220 : ind1 = i - 1;
129 7013220 : ind2 = i;
130 : }
131 :
132 7025300 : double dy1 = padfY[ind1];
133 7025300 : double dy2 = padfY[ind2];
134 :
135 7025300 : if ((dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy))
136 7001160 : continue;
137 :
138 24262 : double dx1 = 0.0;
139 24262 : double dx2 = 0.0;
140 24262 : if (dy1 < dy2)
141 : {
142 12065 : dx1 = padfX[ind1];
143 12065 : dx2 = padfX[ind2];
144 : }
145 12197 : else if (dy1 > dy2)
146 : {
147 12075 : std::swap(dy1, dy2);
148 12075 : dx2 = padfX[ind1];
149 12075 : dx1 = padfX[ind2];
150 : }
151 : else // if( fabs(dy1-dy2) < 1.e-6 )
152 : {
153 :
154 : // AE: DO NOT skip bottom horizontal segments
155 : // -Fill them separately-
156 122 : if (padfX[ind1] > padfX[ind2])
157 : {
158 38 : const double dhorizontal_x1 = floor(padfX[ind2] + 0.5);
159 38 : const double dhorizontal_x2 = floor(padfX[ind1] + 0.5);
160 :
161 38 : if ((dhorizontal_x1 > static_cast<double>(maxx)) ||
162 38 : (dhorizontal_x2 <= 0))
163 0 : continue;
164 : const int horizontal_x1 =
165 38 : static_cast<int>(std::max(dhorizontal_x1, 0.0));
166 : const int horizontal_x2 = static_cast<int>(
167 38 : std::min<double>(dhorizontal_x2, nRasterXSize));
168 :
169 : // Fill the horizontal segment (separately from the rest).
170 38 : if (bAvoidBurningSamePoints)
171 : {
172 10 : polyInts2[ints2++] = horizontal_x1;
173 10 : polyInts2[ints2++] = horizontal_x2;
174 : }
175 : else
176 : {
177 28 : pfnScanlineFunc(
178 : pCBData, y, horizontal_x1, horizontal_x2 - 1,
179 : dfVariant == nullptr ? 0 : dfVariant[0]);
180 : }
181 : }
182 : // else: Skip top horizontal segments.
183 : // They are already filled in the regular loop.
184 122 : continue;
185 : }
186 :
187 24140 : if (dy < dy2 && dy >= dy1)
188 : {
189 : const double intersect = std::clamp<double>(
190 48092 : (dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1, INT_MIN,
191 24046 : INT_MAX);
192 :
193 24046 : polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
194 : }
195 : }
196 :
197 11930 : std::sort(polyInts.begin(), polyInts.begin() + ints);
198 11930 : std::sort(polyInts2.begin(), polyInts2.begin() + ints2);
199 :
200 23953 : for (int i = 0; i + 1 < ints; i += 2)
201 : {
202 12023 : if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
203 : {
204 11312 : pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
205 : dfVariant == nullptr ? 0 : dfVariant[0]);
206 : }
207 : }
208 :
209 11940 : for (int i2 = 0, i = 0; i2 + 1 < ints2; i2 += 2)
210 : {
211 10 : if (polyInts2[i2] <= maxx && polyInts2[i2 + 1] > minx)
212 : {
213 : // "synchronize" polyInts[i] with polyInts2[i2]
214 10 : 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 10 : if (i + 1 >= ints || polyInts[i] != polyInts2[i2])
219 : {
220 4 : pfnScanlineFunc(pCBData, y, polyInts2[i2],
221 2 : polyInts2[i2 + 1] - 1,
222 : dfVariant == nullptr ? 0 : dfVariant[0]);
223 : }
224 : }
225 : }
226 : }
227 : }
228 :
229 : /************************************************************************/
230 : /* GDALdllImagePoint() */
231 : /************************************************************************/
232 :
233 5 : 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 : {
238 10 : for (int i = 0; i < nPartCount; i++)
239 : {
240 5 : const double dfX = padfX[i];
241 5 : const double dfY = padfY[i];
242 5 : double dfVariant = 0.0;
243 5 : if (padfVariant != nullptr)
244 0 : dfVariant = padfVariant[i];
245 :
246 5 : if (0 <= dfX && dfX < nRasterXSize && 0 <= dfY && dfY < nRasterYSize)
247 5 : pfnPointFunc(pCBData, static_cast<int>(dfY), static_cast<int>(dfX),
248 : dfVariant);
249 : }
250 5 : }
251 :
252 : /************************************************************************/
253 : /* GDALdllImageLine() */
254 : /************************************************************************/
255 :
256 1874 : void GDALdllImageLine(int nRasterXSize, int nRasterYSize, int nPartCount,
257 : const int *panPartSize, const double *padfX,
258 : const double *padfY, const double *padfVariant,
259 : llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
260 : {
261 1874 : if (!nPartCount)
262 0 : return;
263 :
264 3748 : for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
265 : {
266 48667 : for (int j = 1; j < panPartSize[i]; j++)
267 : {
268 46793 : double dfX = padfX[n + j - 1];
269 46793 : double dfY = padfY[n + j - 1];
270 46793 : double dfXEnd = padfX[n + j];
271 46793 : double dfYEnd = padfY[n + j];
272 :
273 : // Skip segments that are off the target region.
274 46793 : if ((dfY < 0.0 && dfYEnd < 0.0) ||
275 46793 : (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
276 0 : (dfX < 0.0 && dfXEnd < 0.0) ||
277 46793 : (dfX > nRasterXSize && dfXEnd > nRasterXSize))
278 0 : continue;
279 :
280 : // TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
281 46793 : if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
282 46793 : dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
283 46793 : dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
284 : {
285 0 : CPLErrorOnce(
286 : CE_Warning, CPLE_AppDefined,
287 : "GDALdllImageLine(): coordinates outside of int range");
288 0 : continue;
289 : }
290 :
291 46793 : int iX = static_cast<int>(floor(dfX));
292 46793 : int iY = static_cast<int>(floor(dfY));
293 :
294 46793 : const int iX1 = static_cast<int>(floor(dfXEnd));
295 46793 : const int iY1 = static_cast<int>(floor(dfYEnd));
296 :
297 46793 : double dfVariant = 0.0;
298 46793 : double dfVariant1 = 0.0;
299 46793 : if (padfVariant != nullptr &&
300 46758 : pCBData->eBurnValueSource != GBV_UserBurnValue)
301 : {
302 46758 : dfVariant = padfVariant[n + j - 1];
303 46758 : dfVariant1 = padfVariant[n + j];
304 : }
305 :
306 46793 : int nDeltaX = std::abs(iX1 - iX);
307 46793 : int nDeltaY = std::abs(iY1 - iY);
308 :
309 : // Step direction depends on line direction.
310 46793 : const int nXStep = (iX > iX1) ? -1 : 1;
311 46793 : const int nYStep = (iY > iY1) ? -1 : 1;
312 :
313 : // Determine the line slope.
314 46793 : if (nDeltaX >= nDeltaY)
315 : {
316 38876 : const int nXError = nDeltaY << 1;
317 38876 : const int nYError = nXError - (nDeltaX << 1);
318 38876 : int nError = nXError - nDeltaX;
319 : // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
320 : // but if it is == 0, dfDeltaVariant is not really used, so any
321 : // value is okay.
322 38876 : const double dfDeltaVariant =
323 38876 : nDeltaX == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaX;
324 :
325 : // Do not burn the end point, unless we are in the last
326 : // segment. This is to avoid burning twice intermediate points,
327 : // which causes artifacts in Add mode
328 38876 : if (j != panPartSize[i] - 1)
329 : {
330 37185 : nDeltaX--;
331 : }
332 :
333 67162 : while (nDeltaX-- >= 0)
334 : {
335 28286 : if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
336 : iY < nRasterYSize)
337 28089 : pfnPointFunc(pCBData, iY, iX, dfVariant);
338 :
339 28286 : dfVariant += dfDeltaVariant;
340 28286 : iX += nXStep;
341 28286 : if (nError > 0)
342 : {
343 11024 : iY += nYStep;
344 11024 : nError += nYError;
345 : }
346 : else
347 : {
348 17262 : nError += nXError;
349 : }
350 : }
351 : }
352 : else
353 : {
354 7917 : const int nXError = nDeltaX << 1;
355 7917 : const int nYError = nXError - (nDeltaY << 1);
356 7917 : int nError = nXError - nDeltaY;
357 : // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
358 : // but if it is == 0, dfDeltaVariant is not really used, so any
359 : // value is okay.
360 7917 : double dfDeltaVariant =
361 7917 : nDeltaY == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaY;
362 :
363 : // Do not burn the end point, unless we are in the last
364 : // segment. This is to avoid burning twice intermediate points,
365 : // which causes artifacts in Add mode
366 7917 : if (j != panPartSize[i] - 1)
367 : {
368 7734 : nDeltaY--;
369 : }
370 :
371 16207 : while (nDeltaY-- >= 0)
372 : {
373 8290 : if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
374 : iY < nRasterYSize)
375 8201 : pfnPointFunc(pCBData, iY, iX, dfVariant);
376 :
377 8290 : dfVariant += dfDeltaVariant;
378 8290 : iY += nYStep;
379 8290 : if (nError > 0)
380 : {
381 35 : iX += nXStep;
382 35 : nError += nYError;
383 : }
384 : else
385 : {
386 8255 : nError += nXError;
387 : }
388 : }
389 : }
390 : }
391 : }
392 : }
393 :
394 : /************************************************************************/
395 : /* GDALdllImageLineAllTouched() */
396 : /* */
397 : /* This alternate line drawing algorithm attempts to ensure */
398 : /* that every pixel touched at all by the line will get set. */
399 : /* @param padfVariant should contain the values that are to be */
400 : /* added to the burn value. The values along the line between the */
401 : /* points will be linearly interpolated. These values are used */
402 : /* only if pCBData->eBurnValueSource is set to something other */
403 : /* than GBV_UserBurnValue. If NULL is passed, a monotonous line */
404 : /* will be drawn with the burn value. */
405 : /************************************************************************/
406 :
407 141 : void GDALdllImageLineAllTouched(
408 : int nRasterXSize, int nRasterYSize, int nPartCount, const int *panPartSize,
409 : const double *padfX, const double *padfY, const double *padfVariant,
410 : llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData,
411 : bool bAvoidBurningSamePoints, bool bIntersectOnly)
412 :
413 : {
414 : // This is an epsilon to detect geometries that are aligned with pixel
415 : // coordinates. Hard to find the right value. We put it to that value
416 : // to satisfy the scenarios of https://github.com/OSGeo/gdal/issues/7523
417 : // and https://github.com/OSGeo/gdal/issues/6414
418 141 : constexpr double EPSILON_INTERSECT_ONLY = 1e-4;
419 :
420 141 : if (!nPartCount)
421 0 : return;
422 :
423 285 : for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
424 : {
425 288 : std::set<std::pair<int, int>> lastBurntPoints;
426 288 : std::set<std::pair<int, int>> newBurntPoints;
427 :
428 1000 : for (int j = 1; j < panPartSize[i]; j++)
429 : {
430 856 : lastBurntPoints = std::move(newBurntPoints);
431 856 : newBurntPoints.clear();
432 :
433 856 : double dfX = padfX[n + j - 1];
434 856 : double dfY = padfY[n + j - 1];
435 :
436 856 : double dfXEnd = padfX[n + j];
437 856 : double dfYEnd = padfY[n + j];
438 :
439 : // Skip segments that are off the target region.
440 856 : if ((dfY < 0.0 && dfYEnd < 0.0) ||
441 845 : (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
442 831 : (dfX < 0.0 && dfXEnd < 0.0) ||
443 800 : (dfX > nRasterXSize && dfXEnd > nRasterXSize))
444 553 : continue;
445 :
446 : // TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
447 763 : if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
448 763 : dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
449 763 : dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
450 : {
451 0 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
452 : "GDALdllImageLineAllTouched(): coordinates "
453 : "outside of int range");
454 0 : continue;
455 : }
456 :
457 763 : double dfVariant = 0.0;
458 763 : double dfVariantEnd = 0.0;
459 763 : if (padfVariant != nullptr &&
460 11 : static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
461 : GBV_UserBurnValue)
462 : {
463 11 : dfVariant = padfVariant[n + j - 1];
464 11 : dfVariantEnd = padfVariant[n + j];
465 : }
466 :
467 : // Swap if needed so we can proceed from left2right (X increasing)
468 763 : if (dfX > dfXEnd)
469 : {
470 272 : std::swap(dfX, dfXEnd);
471 272 : std::swap(dfY, dfYEnd);
472 272 : std::swap(dfVariant, dfVariantEnd);
473 : }
474 :
475 : // Special case for vertical lines.
476 :
477 763 : if (fabs(dfX - dfXEnd) < .01)
478 : {
479 209 : if (bIntersectOnly)
480 : {
481 200 : if (std::abs(dfX - std::round(dfX)) <
482 272 : EPSILON_INTERSECT_ONLY &&
483 72 : std::abs(dfXEnd - std::round(dfXEnd)) <
484 : EPSILON_INTERSECT_ONLY)
485 71 : continue;
486 : }
487 :
488 138 : if (dfYEnd < dfY)
489 : {
490 73 : std::swap(dfY, dfYEnd);
491 73 : std::swap(dfVariant, dfVariantEnd);
492 : }
493 :
494 138 : const int iX = static_cast<int>(floor(dfXEnd));
495 138 : int iY = static_cast<int>(floor(dfY));
496 138 : int iYEnd =
497 138 : static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));
498 :
499 138 : if (iX < 0 || iX >= nRasterXSize)
500 1 : continue;
501 :
502 137 : double dfDeltaVariant = 0.0;
503 137 : if (dfYEnd - dfY > 0.0)
504 136 : dfDeltaVariant = (dfVariantEnd - dfVariant) /
505 136 : (dfYEnd - dfY); // Per unit change in iY.
506 :
507 : // Clip to the borders of the target region.
508 137 : if (iY < 0)
509 0 : iY = 0;
510 137 : if (iYEnd >= nRasterYSize)
511 0 : iYEnd = nRasterYSize - 1;
512 137 : dfVariant += dfDeltaVariant * (iY - dfY);
513 :
514 137 : if (padfVariant == nullptr)
515 : {
516 1289 : for (; iY <= iYEnd; iY++)
517 : {
518 1156 : if (bAvoidBurningSamePoints)
519 : {
520 146 : auto yx = std::pair<int, int>(iY, iX);
521 146 : if (lastBurntPoints.find(yx) !=
522 292 : lastBurntPoints.end())
523 : {
524 13 : continue;
525 : }
526 133 : newBurntPoints.insert(yx);
527 : }
528 1143 : pfnPointFunc(pCBData, iY, iX, 0.0);
529 : }
530 : }
531 : else
532 : {
533 14 : for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
534 : {
535 10 : if (bAvoidBurningSamePoints)
536 : {
537 0 : auto yx = std::pair<int, int>(iY, iX);
538 0 : if (lastBurntPoints.find(yx) !=
539 0 : lastBurntPoints.end())
540 : {
541 0 : continue;
542 : }
543 0 : newBurntPoints.insert(yx);
544 : }
545 10 : pfnPointFunc(pCBData, iY, iX, dfVariant);
546 : }
547 : }
548 :
549 137 : continue; // Next segment.
550 : }
551 :
552 554 : const double dfDeltaVariant =
553 554 : (dfVariantEnd - dfVariant) /
554 554 : (dfXEnd - dfX); // Per unit change in iX.
555 :
556 : // Special case for horizontal lines.
557 554 : if (fabs(dfY - dfYEnd) < .01)
558 : {
559 251 : if (bIntersectOnly)
560 : {
561 242 : if (std::abs(dfY - std::round(dfY)) <
562 317 : EPSILON_INTERSECT_ONLY &&
563 75 : std::abs(dfYEnd - std::round(dfYEnd)) <
564 : EPSILON_INTERSECT_ONLY)
565 74 : continue;
566 : }
567 :
568 177 : if (dfXEnd < dfX)
569 : {
570 0 : std::swap(dfX, dfXEnd);
571 0 : std::swap(dfVariant, dfVariantEnd);
572 : }
573 :
574 177 : int iX = static_cast<int>(floor(dfX));
575 177 : const int iY = static_cast<int>(floor(dfY));
576 177 : int iXEnd =
577 177 : static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));
578 :
579 177 : if (iY < 0 || iY >= nRasterYSize)
580 0 : continue;
581 :
582 : // Clip to the borders of the target region.
583 177 : if (iX < 0)
584 2 : iX = 0;
585 177 : if (iXEnd >= nRasterXSize)
586 1 : iXEnd = nRasterXSize - 1;
587 177 : dfVariant += dfDeltaVariant * (iX - dfX);
588 :
589 177 : if (padfVariant == nullptr)
590 : {
591 1163 : for (; iX <= iXEnd; iX++)
592 : {
593 990 : if (bAvoidBurningSamePoints)
594 : {
595 144 : auto yx = std::pair<int, int>(iY, iX);
596 144 : if (lastBurntPoints.find(yx) !=
597 288 : lastBurntPoints.end())
598 : {
599 12 : continue;
600 : }
601 132 : newBurntPoints.insert(yx);
602 : }
603 978 : pfnPointFunc(pCBData, iY, iX, 0.0);
604 : }
605 : }
606 : else
607 : {
608 14 : for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
609 : {
610 10 : if (bAvoidBurningSamePoints)
611 : {
612 0 : auto yx = std::pair<int, int>(iY, iX);
613 0 : if (lastBurntPoints.find(yx) !=
614 0 : lastBurntPoints.end())
615 : {
616 0 : continue;
617 : }
618 0 : newBurntPoints.insert(yx);
619 : }
620 10 : pfnPointFunc(pCBData, iY, iX, dfVariant);
621 : }
622 : }
623 :
624 177 : continue; // Next segment.
625 : }
626 :
627 : /* --------------------------------------------------------------------
628 : */
629 : /* General case - left to right sloped. */
630 : /* --------------------------------------------------------------------
631 : */
632 303 : const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
633 :
634 : // Clip segment in X.
635 303 : if (dfXEnd > nRasterXSize)
636 : {
637 16 : dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
638 16 : dfXEnd = nRasterXSize;
639 : }
640 303 : if (dfX < 0.0)
641 : {
642 18 : dfY += (0.0 - dfX) * dfSlope;
643 18 : dfVariant += dfDeltaVariant * (0.0 - dfX);
644 18 : dfX = 0.0;
645 : }
646 :
647 : // Clip segment in Y.
648 303 : if (dfYEnd > dfY)
649 : {
650 74 : if (dfY < 0.0)
651 : {
652 7 : const double dfDiffX = (0.0 - dfY) / dfSlope;
653 7 : dfX += dfDiffX;
654 7 : dfVariant += dfDeltaVariant * dfDiffX;
655 7 : dfY = 0.0;
656 : }
657 74 : if (dfYEnd >= nRasterYSize)
658 : {
659 9 : dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
660 9 : if (dfXEnd > nRasterXSize)
661 5 : dfXEnd = nRasterXSize;
662 : // dfYEnd is no longer used afterwards, but for
663 : // consistency it should be:
664 : // dfYEnd = nRasterXSize;
665 : }
666 : }
667 : else
668 : {
669 229 : if (dfY >= nRasterYSize)
670 : {
671 13 : const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
672 13 : dfX += dfDiffX;
673 13 : dfVariant += dfDeltaVariant * dfDiffX;
674 13 : dfY = nRasterYSize;
675 : }
676 229 : if (dfYEnd < 0.0)
677 : {
678 7 : dfXEnd -= (dfYEnd - 0) / dfSlope;
679 : // dfYEnd is no longer used afterwards, but for
680 : // consistency it should be:
681 : // dfYEnd = 0.0;
682 : }
683 : }
684 :
685 : // Step from pixel to pixel.
686 6738 : while (dfX >= 0.0 && dfX < dfXEnd)
687 : {
688 6435 : const int iX = static_cast<int>(floor(dfX));
689 6435 : const int iY = static_cast<int>(floor(dfY));
690 :
691 : // Burn in the current point.
692 : // We should be able to drop the Y check because we clipped
693 : // in Y, but there may be some error with all the small steps.
694 6435 : if (iY >= 0 && iY < nRasterYSize)
695 : {
696 6336 : if (bAvoidBurningSamePoints)
697 : {
698 42 : auto yx = std::pair<int, int>(iY, iX);
699 79 : if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
700 79 : newBurntPoints.find(yx) == newBurntPoints.end())
701 : {
702 27 : newBurntPoints.insert(yx);
703 27 : pfnPointFunc(pCBData, iY, iX, dfVariant);
704 : }
705 : }
706 : else
707 : {
708 6294 : pfnPointFunc(pCBData, iY, iX, dfVariant);
709 : }
710 : }
711 :
712 6435 : double dfStepX = floor(dfX + 1.0) - dfX;
713 6435 : double dfStepY = dfStepX * dfSlope;
714 :
715 : // Step to right pixel without changing scanline?
716 6435 : if (static_cast<int>(floor(dfY + dfStepY)) == iY)
717 : {
718 3840 : dfX += dfStepX;
719 3840 : dfY += dfStepY;
720 3840 : dfVariant += dfDeltaVariant * dfStepX;
721 : }
722 2595 : else if (dfSlope < 0)
723 : {
724 1780 : dfStepY = iY - dfY;
725 1780 : if (dfStepY > -0.000000001)
726 877 : dfStepY = -0.000000001;
727 :
728 1780 : dfStepX = dfStepY / dfSlope;
729 1780 : dfX += dfStepX;
730 1780 : dfY += dfStepY;
731 1780 : dfVariant += dfDeltaVariant * dfStepX;
732 : }
733 : else
734 : {
735 815 : dfStepY = (iY + 1) - dfY;
736 815 : if (dfStepY < 0.000000001)
737 1 : dfStepY = 0.000000001;
738 :
739 815 : dfStepX = dfStepY / dfSlope;
740 815 : dfX += dfStepX;
741 815 : dfY += dfStepY;
742 815 : dfVariant += dfDeltaVariant * dfStepX;
743 : }
744 : } // Next step along segment.
745 : } // Next segment.
746 : } // Next part.
747 : }
|