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 490 : static int llCompareInt(const void *a, const void *b)
34 : {
35 490 : return (*(const int *)a) - (*(const int *)b);
36 : }
37 :
38 35 : static void llSwapDouble(double *a, double *b)
39 : {
40 35 : double temp = *a;
41 35 : *a = *b;
42 35 : *b = temp;
43 35 : }
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 20 : 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 20 : if (!nPartCount) {
118 0 : return;
119 : }
120 :
121 20 : n = 0;
122 44 : for( part = 0; part < nPartCount; part++ )
123 24 : n += panPartSize[part];
124 :
125 20 : polyInts = (int *) malloc(sizeof(int) * n);
126 20 : polyAllocated = n;
127 :
128 20 : dminy = padfY[0];
129 20 : dmaxy = padfY[0];
130 115 : for (i=1; (i < n); i++) {
131 :
132 95 : if (padfY[i] < dminy) {
133 23 : dminy = padfY[i];
134 : }
135 95 : if (padfY[i] > dmaxy) {
136 5 : dmaxy = padfY[i];
137 : }
138 : }
139 20 : miny = (int) dminy;
140 20 : maxy = (int) dmaxy;
141 :
142 :
143 20 : if( miny < 0 )
144 0 : miny = 0;
145 20 : if( maxy >= nRasterYSize )
146 1 : maxy = nRasterYSize-1;
147 :
148 :
149 20 : minx = 0;
150 20 : maxx = nRasterXSize - 1;
151 :
152 : /* Fix in 1.3: count a vertex only once */
153 442 : for (y=miny; y <= maxy; y++) {
154 422 : int partoffset = 0;
155 :
156 422 : dy = y +0.5; /* center height of line*/
157 :
158 :
159 422 : part = 0;
160 422 : ints = 0;
161 :
162 : /*Initialize polyInts, otherwise it can sometimes causes a seg fault */
163 422 : memset(polyInts, -1, sizeof(int) * n);
164 :
165 3404 : for (i=0; (i < n); i++) {
166 :
167 :
168 2982 : if( i == partoffset + panPartSize[part] ) {
169 210 : partoffset += panPartSize[part];
170 210 : part++;
171 : }
172 :
173 2982 : if( i == partoffset ) {
174 632 : ind1 = partoffset + panPartSize[part] - 1;
175 632 : ind2 = partoffset;
176 : } else {
177 2350 : ind1 = i-1;
178 2350 : ind2 = i;
179 : }
180 :
181 :
182 2982 : dy1 = padfY[ind1];
183 2982 : dy2 = padfY[ind2];
184 :
185 :
186 2982 : if( (dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy) )
187 2158 : continue;
188 :
189 824 : if (dy1 < dy2) {
190 412 : dx1 = padfX[ind1];
191 412 : dx2 = padfX[ind2];
192 412 : } else if (dy1 > dy2) {
193 412 : dy2 = padfY[ind1];
194 412 : dy1 = padfY[ind2];
195 412 : dx2 = padfX[ind1];
196 412 : 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 824 : if(( dy < dy2 ) && (dy >= dy1))
221 : {
222 :
223 824 : intersect = (dy-dy1) * (dx2-dx1) / (dy2-dy1) + dx1;
224 :
225 824 : 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 422 : qsort(polyInts, ints, sizeof(int), llCompareInt);
237 :
238 :
239 834 : for (i=0; (i < (ints)); i+=2)
240 : {
241 412 : if( polyInts[i] <= maxx && polyInts[i+1] > minx )
242 : {
243 412 : pfnScanlineFunc( pCBData, y, polyInts[i], polyInts[i+1] - 1, (dfVariant == NULL)?0:dfVariant[0] );
244 : }
245 : }
246 : }
247 :
248 20 : free( polyInts );
249 : }
250 :
251 : /************************************************************************/
252 : /* GDALdllImagePoint() */
253 : /************************************************************************/
254 :
255 5 : 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 10 : for ( i = 0; i < nPartCount; i++ )
263 : {
264 5 : int nX = (int)floor( padfX[i] );
265 5 : int nY = (int)floor( padfY[i] );
266 5 : double dfVariant = 0;
267 5 : if( padfVariant != NULL )
268 0 : dfVariant = padfVariant[i];
269 :
270 5 : if ( 0 <= nX && nX < nRasterXSize && 0 <= nY && nY < nRasterYSize )
271 5 : pfnPointFunc( pCBData, nY, nX, dfVariant );
272 : }
273 5 : }
274 :
275 : /************************************************************************/
276 : /* GDALdllImageLine() */
277 : /************************************************************************/
278 :
279 1242 : 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 1242 : if ( !nPartCount )
287 0 : return;
288 :
289 2484 : for ( i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
290 : {
291 : int j;
292 :
293 32420 : for ( j = 1; j < panPartSize[i]; j++ )
294 : {
295 31178 : int iX = (int)floor( padfX[n + j - 1] );
296 31178 : int iY = (int)floor( padfY[n + j - 1] );
297 :
298 31178 : const int iX1 = (int)floor( padfX[n + j] );
299 31178 : const int iY1 = (int)floor( padfY[n + j] );
300 :
301 31178 : double dfVariant = 0, dfVariant1 = 0;
302 31178 : if( padfVariant != NULL &&
303 : ((GDALRasterizeInfo *)pCBData)->eBurnValueSource !=
304 : GBV_UserBurnValue )
305 : {
306 31173 : dfVariant = padfVariant[n + j - 1];
307 31173 : dfVariant1 = padfVariant[n + j];
308 : }
309 :
310 31178 : int nDeltaX = ABS( iX1 - iX );
311 31178 : int nDeltaY = ABS( iY1 - iY );
312 :
313 : // Step direction depends on line direction.
314 31178 : const int nXStep = ( iX > iX1 ) ? -1 : 1;
315 31178 : const int nYStep = ( iY > iY1 ) ? -1 : 1;
316 :
317 : // Determine the line slope.
318 31178 : if ( nDeltaX >= nDeltaY )
319 : {
320 27307 : const int nXError = nDeltaY << 1;
321 27307 : const int nYError = nXError - (nDeltaX << 1);
322 27307 : 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 27307 : (double)nDeltaX;
327 :
328 99035 : while ( nDeltaX-- >= 0 )
329 : {
330 44421 : if ( 0 <= iX && iX < nRasterXSize
331 : && 0 <= iY && iY < nRasterYSize )
332 44246 : pfnPointFunc( pCBData, iY, iX, dfVariant );
333 :
334 44421 : dfVariant += dfDeltaVariant;
335 44421 : iX += nXStep;
336 44421 : if ( nError > 0 )
337 : {
338 15856 : iY += nYStep;
339 15856 : nError += nYError;
340 : }
341 : else
342 28565 : nError += nXError;
343 : }
344 : }
345 : else
346 : {
347 3871 : const int nXError = nDeltaX << 1;
348 3871 : const int nYError = nXError - (nDeltaY << 1);
349 3871 : 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 3871 : (double)nDeltaY;
354 :
355 15542 : while ( nDeltaY-- >= 0 )
356 : {
357 7800 : if ( 0 <= iX && iX < nRasterXSize
358 : && 0 <= iY && iY < nRasterYSize )
359 7712 : pfnPointFunc( pCBData, iY, iX, dfVariant );
360 :
361 7800 : dfVariant += dfDeltaVariant;
362 7800 : iY += nYStep;
363 7800 : if ( nError > 0 )
364 : {
365 30 : iX += nXStep;
366 30 : nError += nYError;
367 : }
368 : else
369 7770 : 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 7 : 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 7 : if ( !nPartCount )
399 0 : return;
400 :
401 14 : for ( i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
402 : {
403 : int j;
404 :
405 33 : for ( j = 1; j < panPartSize[i]; j++ )
406 : {
407 26 : double dfX = padfX[n + j - 1];
408 26 : double dfY = padfY[n + j - 1];
409 :
410 26 : double dfXEnd = padfX[n + j];
411 26 : double dfYEnd = padfY[n + j];
412 :
413 26 : double dfVariant = 0, dfVariantEnd = 0;
414 26 : 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 26 : 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 26 : if( dfX > dfXEnd )
431 : {
432 7 : llSwapDouble( &dfX, &dfXEnd );
433 7 : llSwapDouble( &dfY, &dfYEnd );
434 7 : llSwapDouble( &dfVariant, &dfVariantEnd );
435 : }
436 :
437 : // Special case for vertical lines.
438 26 : if( floor(dfX) == floor(dfXEnd) )
439 : {
440 12 : if( dfYEnd < dfY )
441 : {
442 7 : llSwapDouble( &dfY, &dfYEnd );
443 7 : llSwapDouble( &dfVariant, &dfVariantEnd );
444 : }
445 :
446 12 : int iX = (int) floor(dfX);
447 12 : int iY = (int) floor(dfY);
448 12 : int iYEnd = (int) floor(dfYEnd);
449 :
450 12 : if( iX >= nRasterXSize )
451 0 : continue;
452 :
453 12 : double dfDeltaVariant = 0;
454 12 : if(( dfYEnd - dfY ) > 0)
455 : dfDeltaVariant = ( dfVariantEnd - dfVariant )
456 12 : / ( dfYEnd - dfY );//per unit change in iY
457 :
458 : // Clip to the borders of the target region
459 12 : if( iY < 0 )
460 0 : iY = 0;
461 12 : if( iYEnd >= nRasterYSize )
462 0 : iYEnd = nRasterYSize - 1;
463 12 : dfVariant += dfDeltaVariant * ( (double)iY - dfY );
464 :
465 12 : if( padfVariant == NULL )
466 112 : for( ; iY <= iYEnd; iY++ )
467 100 : 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 12 : continue; // next segment
473 : }
474 :
475 : double dfDeltaVariant = ( dfVariantEnd - dfVariant )
476 14 : / ( dfXEnd - dfX );//per unit change in iX
477 :
478 : // special case for horizontal lines
479 14 : if( floor(dfY) == floor(dfYEnd) )
480 : {
481 12 : if( dfXEnd < dfX )
482 : {
483 0 : llSwapDouble( &dfX, &dfXEnd );
484 0 : llSwapDouble( &dfVariant, &dfVariantEnd );
485 : }
486 :
487 12 : int iX = (int) floor(dfX);
488 12 : int iY = (int) floor(dfY);
489 12 : int iXEnd = (int) floor(dfXEnd);
490 :
491 12 : if( iY >= nRasterYSize )
492 0 : continue;
493 :
494 : // Clip to the borders of the target region
495 12 : if( iX < 0 )
496 0 : iX = 0;
497 12 : if( iXEnd >= nRasterXSize )
498 0 : iXEnd = nRasterXSize - 1;
499 12 : dfVariant += dfDeltaVariant * ( (double)iX - dfX );
500 :
501 12 : if( padfVariant == NULL )
502 136 : for( ; iX <= iXEnd; iX++ )
503 124 : 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 12 : continue; // next segment
509 : }
510 :
511 : /* -------------------------------------------------------------------- */
512 : /* General case - left to right sloped. */
513 : /* -------------------------------------------------------------------- */
514 2 : double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
515 :
516 : // clip segment in X
517 2 : if( dfXEnd > nRasterXSize )
518 : {
519 0 : dfYEnd -= ( dfXEnd - (double)nRasterXSize ) * dfSlope;
520 0 : dfXEnd = nRasterXSize;
521 : }
522 2 : 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 2 : 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 2 : if( dfY >= nRasterYSize )
548 : {
549 0 : dfX += (dfDiffX = ((double)nRasterYSize - dfY) / dfSlope);
550 0 : dfVariant += dfDeltaVariant * dfDiffX;
551 0 : dfY = nRasterYSize;
552 : }
553 2 : if( dfYEnd < 0 )
554 : {
555 0 : dfXEnd -= ( dfYEnd - 0 ) / dfSlope;
556 0 : dfYEnd = 0;
557 : }
558 : }
559 :
560 : // step from pixel to pixel.
561 14 : while( dfX < dfXEnd )
562 : {
563 10 : int iX = (int) floor(dfX);
564 10 : 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 10 : if( iY >= 0 && iY < nRasterYSize )
570 10 : pfnPointFunc( pCBData, iY, iX, dfVariant );
571 :
572 10 : double dfStepX = floor(dfX+1.0) - dfX;
573 10 : double dfStepY = dfStepX * dfSlope;
574 :
575 : // step to right pixel without changing scanline?
576 10 : if( (int) floor(dfY + dfStepY) == iY )
577 : {
578 6 : dfX += dfStepX;
579 6 : dfY += dfStepY;
580 6 : dfVariant += dfDeltaVariant * dfStepX;
581 : }
582 4 : else if( dfSlope < 0 )
583 : {
584 4 : dfStepY = iY - dfY;
585 4 : if( dfStepY > -0.000000001 )
586 4 : dfStepY = -0.000000001;
587 :
588 4 : dfStepX = dfStepY / dfSlope;
589 4 : dfX += dfStepX;
590 4 : dfY += dfStepY;
591 4 : 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 :
|
