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