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