1 : /******************************************************************************
2 : * $Id: llrasterize.cpp 22998 2011-08-28 12:23:29Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Vector polygon rasterization code.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.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 "gdal_alg.h"
31 : #include "gdal_alg_priv.h"
32 :
33 980 : static int llCompareInt(const void *a, const void *b)
34 : {
35 980 : return (*(const int *)a) - (*(const int *)b);
36 : }
37 :
38 70 : static void llSwapDouble(double *a, double *b)
39 : {
40 70 : double temp = *a;
41 70 : *a = *b;
42 70 : *b = temp;
43 70 : }
44 :
45 : /************************************************************************/
46 : /* dllImageFilledPolygon() */
47 : /* */
48 : /* Perform scanline conversion of the passed multi-ring */
49 : /* polygon. Note the polygon does not need to be explicitly */
50 : /* closed. The scanline function will be called with */
51 : /* horizontal scanline chunks which may not be entirely */
52 : /* contained within the valid raster area (in the X */
53 : /* direction). */
54 : /* */
55 : /* NEW: Nodes' coordinate are kept as double in order */
56 : /* to compute accurately the intersections with the lines */
57 : /* */
58 : /* Two methods for determining the border pixels: */
59 : /* */
60 : /* 1) method = 0 */
61 : /* Inherits algorithm from version above but with several bugs */
62 : /* fixed except for the cone facing down. */
63 : /* A pixel on which a line intersects a segment of a */
64 : /* polygon will always be considered as inside the shape. */
65 : /* Note that we only compute intersections with lines that */
66 : /* passes through the middle of a pixel (line coord = 0.5, */
67 : /* 1.5, 2.5, etc.) */
68 : /* */
69 : /* 2) method = 1: */
70 : /* A pixel is considered inside a polygon if its center */
71 : /* falls inside the polygon. This is somehow more robust unless */
72 : /* the nodes are placed in the center of the pixels in which */
73 : /* case, due to numerical inaccuracies, it's hard to predict */
74 : /* if the pixel will be considered inside or outside the shape. */
75 : /************************************************************************/
76 :
77 : /*
78 : * NOTE: This code was originally adapted from the gdImageFilledPolygon()
79 : * function in libgd.
80 : *
81 : * http://www.boutell.com/gd/
82 : *
83 : * It was later adapted for direct inclusion in GDAL and relicensed under
84 : * the GDAL MIT/X license (pulled from the OpenEV distribution).
85 : */
86 :
87 40 : void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
88 : int nPartCount, int *panPartSize,
89 : double *padfX, double *padfY,
90 : double *dfVariant,
91 : llScanlineFunc pfnScanlineFunc, void *pCBData )
92 : {
93 : /*************************************************************************
94 : 2nd Method (method=1):
95 : =====================
96 : No known bug
97 : *************************************************************************/
98 :
99 : int i;
100 : int y;
101 : int miny, maxy,minx,maxx;
102 : double dminy, dmaxy;
103 : double dx1, dy1;
104 : double dx2, dy2;
105 : double dy;
106 : double intersect;
107 :
108 :
109 : int ind1, ind2;
110 : int ints, n, part;
111 : int *polyInts, polyAllocated;
112 :
113 :
114 : int horizontal_x1, horizontal_x2;
115 :
116 :
117 40 : if (!nPartCount) {
118 0 : return;
119 : }
120 :
121 40 : n = 0;
122 88 : for( part = 0; part < nPartCount; part++ )
123 48 : n += panPartSize[part];
124 :
125 40 : polyInts = (int *) malloc(sizeof(int) * n);
126 40 : polyAllocated = n;
127 :
128 40 : dminy = padfY[0];
129 40 : dmaxy = padfY[0];
130 230 : for (i=1; (i < n); i++) {
131 :
132 190 : if (padfY[i] < dminy) {
133 46 : dminy = padfY[i];
134 : }
135 190 : if (padfY[i] > dmaxy) {
136 10 : dmaxy = padfY[i];
137 : }
138 : }
139 40 : miny = (int) dminy;
140 40 : maxy = (int) dmaxy;
141 :
142 :
143 40 : if( miny < 0 )
144 0 : miny = 0;
145 40 : if( maxy >= nRasterYSize )
146 2 : maxy = nRasterYSize-1;
147 :
148 :
149 40 : minx = 0;
150 40 : maxx = nRasterXSize - 1;
151 :
152 : /* Fix in 1.3: count a vertex only once */
153 884 : for (y=miny; y <= maxy; y++) {
154 844 : int partoffset = 0;
155 :
156 844 : dy = y +0.5; /* center height of line*/
157 :
158 :
159 844 : part = 0;
160 844 : ints = 0;
161 :
162 : /*Initialize polyInts, otherwise it can sometimes causes a seg fault */
163 844 : memset(polyInts, -1, sizeof(int) * n);
164 :
165 6808 : for (i=0; (i < n); i++) {
166 :
167 :
168 5964 : if( i == partoffset + panPartSize[part] ) {
169 420 : partoffset += panPartSize[part];
170 420 : part++;
171 : }
172 :
173 5964 : if( i == partoffset ) {
174 1264 : ind1 = partoffset + panPartSize[part] - 1;
175 1264 : ind2 = partoffset;
176 : } else {
177 4700 : ind1 = i-1;
178 4700 : ind2 = i;
179 : }
180 :
181 :
182 5964 : dy1 = padfY[ind1];
183 5964 : dy2 = padfY[ind2];
184 :
185 :
186 5964 : if( (dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy) )
187 4316 : continue;
188 :
189 1648 : if (dy1 < dy2) {
190 824 : dx1 = padfX[ind1];
191 824 : dx2 = padfX[ind2];
192 824 : } else if (dy1 > dy2) {
193 824 : dy2 = padfY[ind1];
194 824 : dy1 = padfY[ind2];
195 824 : dx2 = padfX[ind1];
196 824 : dx1 = padfX[ind2];
197 : } else /* if (fabs(dy1-dy2)< 1.e-6) */
198 : {
199 :
200 : /*AE: DO NOT skip bottom horizontal segments
201 : -Fill them separately-
202 : They are not taken into account twice.*/
203 0 : if (padfX[ind1] > padfX[ind2])
204 : {
205 0 : horizontal_x1 = (int) floor(padfX[ind2]+0.5);
206 0 : horizontal_x2 = (int) floor(padfX[ind1]+0.5);
207 :
208 0 : if ( (horizontal_x1 > maxx) || (horizontal_x2 <= minx) )
209 0 : continue;
210 :
211 : /*fill the horizontal segment (separately from the rest)*/
212 0 : pfnScanlineFunc( pCBData, y, horizontal_x1, horizontal_x2 - 1, (dfVariant == NULL)?0:dfVariant[0] );
213 0 : continue;
214 : }
215 : else /*skip top horizontal segments (they are already filled in the regular loop)*/
216 0 : continue;
217 :
218 : }
219 :
220 1648 : if(( dy < dy2 ) && (dy >= dy1))
221 : {
222 :
223 1648 : intersect = (dy-dy1) * (dx2-dx1) / (dy2-dy1) + dx1;
224 :
225 1648 : polyInts[ints++] = (int) floor(intersect+0.5);
226 : }
227 : }
228 :
229 : /*
230 : * It would be more efficient to do this inline, to avoid
231 : * a function call for each comparison.
232 : * NOTE - mloskot: make llCompareInt a functor and use std
233 : * algorithm and it will be optimized and expanded
234 : * automatically in compile-time, with modularity preserved.
235 : */
236 844 : qsort(polyInts, ints, sizeof(int), llCompareInt);
237 :
238 :
239 1668 : for (i=0; (i < (ints)); i+=2)
240 : {
241 824 : if( polyInts[i] <= maxx && polyInts[i+1] > minx )
242 : {
243 824 : pfnScanlineFunc( pCBData, y, polyInts[i], polyInts[i+1] - 1, (dfVariant == NULL)?0:dfVariant[0] );
244 : }
245 : }
246 : }
247 :
248 40 : free( polyInts );
249 : }
250 :
251 : /************************************************************************/
252 : /* GDALdllImagePoint() */
253 : /************************************************************************/
254 :
255 10 : void GDALdllImagePoint( int nRasterXSize, int nRasterYSize,
256 : int nPartCount, int *panPartSize,
257 : double *padfX, double *padfY, double *padfVariant,
258 : llPointFunc pfnPointFunc, void *pCBData )
259 : {
260 : int i;
261 :
262 20 : for ( i = 0; i < nPartCount; i++ )
263 : {
264 10 : int nX = (int)floor( padfX[i] );
265 10 : int nY = (int)floor( padfY[i] );
266 10 : double dfVariant = 0;
267 10 : if( padfVariant != NULL )
268 0 : dfVariant = padfVariant[i];
269 :
270 10 : if ( 0 <= nX && nX < nRasterXSize && 0 <= nY && nY < nRasterYSize )
271 10 : pfnPointFunc( pCBData, nY, nX, dfVariant );
272 : }
273 10 : }
274 :
275 : /************************************************************************/
276 : /* GDALdllImageLine() */
277 : /************************************************************************/
278 :
279 2484 : void GDALdllImageLine( int nRasterXSize, int nRasterYSize,
280 : int nPartCount, int *panPartSize,
281 : double *padfX, double *padfY, double *padfVariant,
282 : llPointFunc pfnPointFunc, void *pCBData )
283 : {
284 : int i, n;
285 :
286 2484 : if ( !nPartCount )
287 0 : return;
288 :
289 4968 : for ( i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
290 : {
291 : int j;
292 :
293 64840 : for ( j = 1; j < panPartSize[i]; j++ )
294 : {
295 62356 : int iX = (int)floor( padfX[n + j - 1] );
296 62356 : int iY = (int)floor( padfY[n + j - 1] );
297 :
298 62356 : const int iX1 = (int)floor( padfX[n + j] );
299 62356 : const int iY1 = (int)floor( padfY[n + j] );
300 :
301 62356 : double dfVariant = 0, dfVariant1 = 0;
302 62356 : if( padfVariant != NULL &&
303 : ((GDALRasterizeInfo *)pCBData)->eBurnValueSource !=
304 : GBV_UserBurnValue )
305 : {
306 62346 : dfVariant = padfVariant[n + j - 1];
307 62346 : dfVariant1 = padfVariant[n + j];
308 : }
309 :
310 62356 : int nDeltaX = ABS( iX1 - iX );
311 62356 : int nDeltaY = ABS( iY1 - iY );
312 :
313 : // Step direction depends on line direction.
314 62356 : const int nXStep = ( iX > iX1 ) ? -1 : 1;
315 62356 : const int nYStep = ( iY > iY1 ) ? -1 : 1;
316 :
317 : // Determine the line slope.
318 62356 : if ( nDeltaX >= nDeltaY )
319 : {
320 54614 : const int nXError = nDeltaY << 1;
321 54614 : const int nYError = nXError - (nDeltaX << 1);
322 54614 : int nError = nXError - nDeltaX;
323 : /* == 0 makes clang -fcatch-undefined-behavior -ftrapv happy, but if */
324 : /* it is == 0, dfDeltaVariant is not really used, so any value is OK */
325 : double dfDeltaVariant = (nDeltaX == 0) ? 0 : (dfVariant1 - dfVariant) /
326 54614 : (double)nDeltaX;
327 :
328 198070 : while ( nDeltaX-- >= 0 )
329 : {
330 88842 : if ( 0 <= iX && iX < nRasterXSize
331 : && 0 <= iY && iY < nRasterYSize )
332 88492 : pfnPointFunc( pCBData, iY, iX, dfVariant );
333 :
334 88842 : dfVariant += dfDeltaVariant;
335 88842 : iX += nXStep;
336 88842 : if ( nError > 0 )
337 : {
338 31712 : iY += nYStep;
339 31712 : nError += nYError;
340 : }
341 : else
342 57130 : nError += nXError;
343 : }
344 : }
345 : else
346 : {
347 7742 : const int nXError = nDeltaX << 1;
348 7742 : const int nYError = nXError - (nDeltaY << 1);
349 7742 : int nError = nXError - nDeltaY;
350 : /* == 0 makes clang -fcatch-undefined-behavior -ftrapv happy, but if */
351 : /* it is == 0, dfDeltaVariant is not really used, so any value is OK */
352 : double dfDeltaVariant = (nDeltaY == 0) ? 0 : (dfVariant1 - dfVariant) /
353 7742 : (double)nDeltaY;
354 :
355 31084 : while ( nDeltaY-- >= 0 )
356 : {
357 15600 : if ( 0 <= iX && iX < nRasterXSize
358 : && 0 <= iY && iY < nRasterYSize )
359 15424 : pfnPointFunc( pCBData, iY, iX, dfVariant );
360 :
361 15600 : dfVariant += dfDeltaVariant;
362 15600 : iY += nYStep;
363 15600 : if ( nError > 0 )
364 : {
365 60 : iX += nXStep;
366 60 : nError += nYError;
367 : }
368 : else
369 15540 : nError += nXError;
370 : }
371 : }
372 : }
373 : }
374 : }
375 :
376 : /************************************************************************/
377 : /* GDALdllImageLineAllTouched() */
378 : /* */
379 : /* This alternate line drawing algorithm attempts to ensure */
380 : /* that every pixel touched at all by the line will get set. */
381 : /* @param padfVariant should contain the values that are to be */
382 : /* added to the burn value. The values along the line between the */
383 : /* points will be linearly interpolated. These values are used */
384 : /* only if pCBData->eBurnValueSource is set to something other */
385 : /* than GBV_UserBurnValue. If NULL is passed, a monotonous line */
386 : /* will be drawn with the burn value. */
387 : /************************************************************************/
388 :
389 : void
390 14 : GDALdllImageLineAllTouched(int nRasterXSize, int nRasterYSize,
391 : int nPartCount, int *panPartSize,
392 : double *padfX, double *padfY, double *padfVariant,
393 : llPointFunc pfnPointFunc, void *pCBData )
394 :
395 : {
396 : int i, n;
397 :
398 14 : if ( !nPartCount )
399 0 : return;
400 :
401 28 : for ( i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
402 : {
403 : int j;
404 :
405 66 : for ( j = 1; j < panPartSize[i]; j++ )
406 : {
407 52 : double dfX = padfX[n + j - 1];
408 52 : double dfY = padfY[n + j - 1];
409 :
410 52 : double dfXEnd = padfX[n + j];
411 52 : double dfYEnd = padfY[n + j];
412 :
413 52 : double dfVariant = 0, dfVariantEnd = 0;
414 52 : if( padfVariant != NULL &&
415 : ((GDALRasterizeInfo *)pCBData)->eBurnValueSource !=
416 : GBV_UserBurnValue )
417 : {
418 0 : dfVariant = padfVariant[n + j - 1];
419 0 : dfVariantEnd = padfVariant[n + j];
420 : }
421 :
422 : // Skip segments that are off the target region.
423 52 : if( (dfY < 0 && dfYEnd < 0)
424 : || (dfY > nRasterYSize && dfYEnd > nRasterYSize)
425 : || (dfX < 0 && dfXEnd < 0)
426 : || (dfX > nRasterXSize && dfXEnd > nRasterXSize) )
427 0 : continue;
428 :
429 : // Swap if needed so we can proceed from left2right (X increasing)
430 52 : if( dfX > dfXEnd )
431 : {
432 14 : llSwapDouble( &dfX, &dfXEnd );
433 14 : llSwapDouble( &dfY, &dfYEnd );
434 14 : llSwapDouble( &dfVariant, &dfVariantEnd );
435 : }
436 :
437 : // Special case for vertical lines.
438 52 : if( floor(dfX) == floor(dfXEnd) )
439 : {
440 24 : if( dfYEnd < dfY )
441 : {
442 14 : llSwapDouble( &dfY, &dfYEnd );
443 14 : llSwapDouble( &dfVariant, &dfVariantEnd );
444 : }
445 :
446 24 : int iX = (int) floor(dfX);
447 24 : int iY = (int) floor(dfY);
448 24 : int iYEnd = (int) floor(dfYEnd);
449 :
450 24 : if( iX >= nRasterXSize )
451 0 : continue;
452 :
453 24 : double dfDeltaVariant = 0;
454 24 : if(( dfYEnd - dfY ) > 0)
455 : dfDeltaVariant = ( dfVariantEnd - dfVariant )
456 24 : / ( dfYEnd - dfY );//per unit change in iY
457 :
458 : // Clip to the borders of the target region
459 24 : if( iY < 0 )
460 0 : iY = 0;
461 24 : if( iYEnd >= nRasterYSize )
462 0 : iYEnd = nRasterYSize - 1;
463 24 : dfVariant += dfDeltaVariant * ( (double)iY - dfY );
464 :
465 24 : if( padfVariant == NULL )
466 224 : for( ; iY <= iYEnd; iY++ )
467 200 : pfnPointFunc( pCBData, iY, iX, 0.0 );
468 : else
469 0 : for( ; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant )
470 0 : pfnPointFunc( pCBData, iY, iX, dfVariant );
471 :
472 24 : continue; // next segment
473 : }
474 :
475 : double dfDeltaVariant = ( dfVariantEnd - dfVariant )
476 28 : / ( dfXEnd - dfX );//per unit change in iX
477 :
478 : // special case for horizontal lines
479 28 : if( floor(dfY) == floor(dfYEnd) )
480 : {
481 24 : if( dfXEnd < dfX )
482 : {
483 0 : llSwapDouble( &dfX, &dfXEnd );
484 0 : llSwapDouble( &dfVariant, &dfVariantEnd );
485 : }
486 :
487 24 : int iX = (int) floor(dfX);
488 24 : int iY = (int) floor(dfY);
489 24 : int iXEnd = (int) floor(dfXEnd);
490 :
491 24 : if( iY >= nRasterYSize )
492 0 : continue;
493 :
494 : // Clip to the borders of the target region
495 24 : if( iX < 0 )
496 0 : iX = 0;
497 24 : if( iXEnd >= nRasterXSize )
498 0 : iXEnd = nRasterXSize - 1;
499 24 : dfVariant += dfDeltaVariant * ( (double)iX - dfX );
500 :
501 24 : if( padfVariant == NULL )
502 272 : for( ; iX <= iXEnd; iX++ )
503 248 : pfnPointFunc( pCBData, iY, iX, 0.0 );
504 : else
505 0 : for( ; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant )
506 0 : pfnPointFunc( pCBData, iY, iX, dfVariant );
507 :
508 24 : continue; // next segment
509 : }
510 :
511 : /* -------------------------------------------------------------------- */
512 : /* General case - left to right sloped. */
513 : /* -------------------------------------------------------------------- */
514 4 : double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
515 :
516 : // clip segment in X
517 4 : if( dfXEnd > nRasterXSize )
518 : {
519 0 : dfYEnd -= ( dfXEnd - (double)nRasterXSize ) * dfSlope;
520 0 : dfXEnd = nRasterXSize;
521 : }
522 4 : if( dfX < 0 )
523 : {
524 0 : dfY += (0 - dfX) * dfSlope;
525 0 : dfVariant += dfDeltaVariant * (0.0 - dfX);
526 0 : dfX = 0.0;
527 : }
528 :
529 : // clip segment in Y
530 : double dfDiffX;
531 4 : if( dfYEnd > dfY )
532 : {
533 0 : if( dfY < 0 )
534 : {
535 0 : dfX += (dfDiffX = (0 - dfY) / dfSlope);
536 0 : dfVariant += dfDeltaVariant * dfDiffX;
537 0 : dfY = 0.0;
538 : }
539 0 : if( dfYEnd >= nRasterYSize )
540 : {
541 0 : dfXEnd += ( dfYEnd - (double)nRasterYSize ) / dfSlope;
542 0 : dfYEnd = nRasterXSize;
543 : }
544 : }
545 : else
546 : {
547 4 : if( dfY >= nRasterYSize )
548 : {
549 0 : dfX += (dfDiffX = ((double)nRasterYSize - dfY) / dfSlope);
550 0 : dfVariant += dfDeltaVariant * dfDiffX;
551 0 : dfY = nRasterYSize;
552 : }
553 4 : if( dfYEnd < 0 )
554 : {
555 0 : dfXEnd -= ( dfYEnd - 0 ) / dfSlope;
556 0 : dfYEnd = 0;
557 : }
558 : }
559 :
560 : // step from pixel to pixel.
561 28 : while( dfX < dfXEnd )
562 : {
563 20 : int iX = (int) floor(dfX);
564 20 : int iY = (int) floor(dfY);
565 :
566 : // burn in the current point.
567 : // We should be able to drop the Y check because we cliped in Y,
568 : // but there may be some error with all the small steps.
569 20 : if( iY >= 0 && iY < nRasterYSize )
570 20 : pfnPointFunc( pCBData, iY, iX, dfVariant );
571 :
572 20 : double dfStepX = floor(dfX+1.0) - dfX;
573 20 : double dfStepY = dfStepX * dfSlope;
574 :
575 : // step to right pixel without changing scanline?
576 20 : if( (int) floor(dfY + dfStepY) == iY )
577 : {
578 12 : dfX += dfStepX;
579 12 : dfY += dfStepY;
580 12 : dfVariant += dfDeltaVariant * dfStepX;
581 : }
582 8 : else if( dfSlope < 0 )
583 : {
584 8 : dfStepY = iY - dfY;
585 8 : if( dfStepY > -0.000000001 )
586 8 : dfStepY = -0.000000001;
587 :
588 8 : dfStepX = dfStepY / dfSlope;
589 8 : dfX += dfStepX;
590 8 : dfY += dfStepY;
591 8 : dfVariant += dfDeltaVariant * dfStepX;
592 : }
593 : else
594 : {
595 0 : dfStepY = (iY + 1) - dfY;
596 0 : if( dfStepY < 0.000000001 )
597 0 : dfStepY = 0.000000001;
598 :
599 0 : dfStepX = dfStepY / dfSlope;
600 0 : dfX += dfStepX;
601 0 : dfY += dfStepY;
602 0 : dfVariant += dfDeltaVariant * dfStepX;
603 : }
604 : } // next step along sement.
605 :
606 : } // next segment
607 : } // next part
608 : }
609 :
|
