1 : /******************************************************************************
2 : * $Id: llrasterize.cpp 18340 2009-12-19 04:41:45Z chaitanya $
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 415 : static int llCompareInt(const void *a, const void *b)
34 : {
35 415 : return (*(const int *)a) - (*(const int *)b);
36 : }
37 :
38 41 : static void llSwapDouble(double *a, double *b)
39 : {
40 41 : double temp = *a;
41 41 : *a = *b;
42 41 : *b = temp;
43 41 : }
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 17 : 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 17 : if (!nPartCount) {
118 0 : return;
119 : }
120 :
121 17 : n = 0;
122 37 : for( part = 0; part < nPartCount; part++ )
123 20 : n += panPartSize[part];
124 :
125 17 : polyInts = (int *) malloc(sizeof(int) * n);
126 17 : polyAllocated = n;
127 :
128 17 : dminy = padfY[0];
129 17 : dmaxy = padfY[0];
130 95 : for (i=1; (i < n); i++) {
131 :
132 78 : if (padfY[i] < dminy) {
133 18 : dminy = padfY[i];
134 : }
135 78 : if (padfY[i] > dmaxy) {
136 8 : dmaxy = padfY[i];
137 : }
138 : }
139 17 : miny = (int) dminy;
140 17 : maxy = (int) dmaxy;
141 :
142 :
143 17 : if( miny < 0 )
144 0 : miny = 0;
145 17 : if( maxy >= nRasterYSize )
146 1 : maxy = nRasterYSize-1;
147 :
148 :
149 17 : minx = 0;
150 17 : maxx = nRasterXSize - 1;
151 :
152 : /* Fix in 1.3: count a vertex only once */
153 393 : for (y=miny; y <= maxy; y++) {
154 376 : int partoffset = 0;
155 :
156 376 : dy = y +0.5; /* center height of line*/
157 :
158 :
159 376 : part = 0;
160 376 : ints = 0;
161 :
162 : /*Initialize polyInts, otherwise it can sometimes causes a seg fault */
163 376 : memset(polyInts, -1, sizeof(int) * n);
164 :
165 3053 : for (i=0; (i < n); i++) {
166 :
167 :
168 2677 : if( i == partoffset + panPartSize[part] ) {
169 195 : partoffset += panPartSize[part];
170 195 : part++;
171 : }
172 :
173 2677 : if( i == partoffset ) {
174 571 : ind1 = partoffset + panPartSize[part] - 1;
175 571 : ind2 = partoffset;
176 : } else {
177 2106 : ind1 = i-1;
178 2106 : ind2 = i;
179 : }
180 :
181 :
182 2677 : dy1 = padfY[ind1];
183 2677 : dy2 = padfY[ind2];
184 :
185 :
186 2677 : if( (dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy) )
187 1955 : continue;
188 :
189 722 : if (dy1 < dy2) {
190 361 : dx1 = padfX[ind1];
191 361 : dx2 = padfX[ind2];
192 361 : } else if (dy1 > dy2) {
193 361 : dy2 = padfY[ind1];
194 361 : dy1 = padfY[ind2];
195 361 : dx2 = padfX[ind1];
196 361 : 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 722 : if(( dy < dy2 ) && (dy >= dy1))
221 : {
222 :
223 722 : intersect = (dy-dy1) * (dx2-dx1) / (dy2-dy1) + dx1;
224 :
225 722 : 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 376 : qsort(polyInts, ints, sizeof(int), llCompareInt);
237 :
238 :
239 737 : for (i=0; (i < (ints)); i+=2)
240 : {
241 361 : if( polyInts[i] <= maxx && polyInts[i+1] > minx )
242 : {
243 361 : pfnScanlineFunc( pCBData, y, polyInts[i], polyInts[i+1] - 1, (dfVariant == NULL)?0:dfVariant[0] );
244 : }
245 : }
246 : }
247 :
248 17 : free( polyInts );
249 : }
250 :
251 : /************************************************************************/
252 : /* GDALdllImagePoint() */
253 : /************************************************************************/
254 :
255 0 : 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 0 : for ( i = 0; i < nPartCount; i++ )
263 : {
264 0 : int nX = (int)floor( padfX[i] + 0.5 );
265 0 : int nY = (int)floor( padfY[i] + 0.5 );
266 0 : double dfVariant = 0;
267 0 : if( padfVariant != NULL )
268 0 : dfVariant = padfVariant[i];
269 :
270 0 : if ( 0 <= nX && nX < nRasterXSize && 0 <= nY && nY < nRasterYSize )
271 0 : pfnPointFunc( pCBData, nY, nX, dfVariant );
272 : }
273 0 : }
274 :
275 : /************************************************************************/
276 : /* GDALdllImageLine() */
277 : /************************************************************************/
278 :
279 4 : 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 4 : if ( !nPartCount )
287 0 : return;
288 :
289 8 : for ( i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
290 : {
291 : int j;
292 :
293 9 : for ( j = 1; j < panPartSize[i]; j++ )
294 : {
295 5 : int iX = (int)floor( padfX[n + j - 1] + 0.5 );
296 5 : int iY = (int)floor( padfY[n + j - 1] + 0.5 );
297 :
298 5 : const int iX1 = (int)floor( padfX[n + j] + 0.5 );
299 5 : const int iY1 = (int)floor( padfY[n + j] + 0.5 );
300 :
301 5 : double dfVariant = 0, dfVariant1 = 0;
302 5 : if( padfVariant != NULL &&
303 : ((GDALRasterizeInfo *)pCBData)->eBurnValueSource !=
304 : GBV_UserBurnValue )
305 : {
306 3 : dfVariant = padfVariant[n + j - 1];
307 3 : dfVariant1 = padfVariant[n + j];
308 : }
309 :
310 5 : int nDeltaX = ABS( iX1 - iX );
311 5 : int nDeltaY = ABS( iY1 - iY );
312 :
313 : // Step direction depends on line direction.
314 5 : const int nXStep = ( iX > iX1 ) ? -1 : 1;
315 5 : const int nYStep = ( iY > iY1 ) ? -1 : 1;
316 :
317 : // Determine the line slope.
318 5 : if ( nDeltaX >= nDeltaY )
319 : {
320 4 : const int nXError = nDeltaY << 1;
321 4 : const int nYError = nXError - (nDeltaX << 1);
322 4 : int nError = nXError - nDeltaX;
323 : double dfDeltaVariant = (dfVariant1 - dfVariant) /
324 4 : (double)nDeltaX;
325 :
326 402 : while ( nDeltaX-- >= 0 )
327 : {
328 394 : if ( 0 <= iX && iX < nRasterXSize
329 : && 0 <= iY && iY < nRasterYSize )
330 387 : pfnPointFunc( pCBData, iY, iX, dfVariant );
331 :
332 394 : dfVariant += dfDeltaVariant;
333 394 : iX += nXStep;
334 394 : if ( nError > 0 )
335 : {
336 202 : iY += nYStep;
337 202 : nError += nYError;
338 : }
339 : else
340 192 : nError += nXError;
341 : }
342 : }
343 : else
344 : {
345 1 : const int nXError = nDeltaX << 1;
346 1 : const int nYError = nXError - (nDeltaY << 1);
347 1 : int nError = nXError - nDeltaY;
348 : double dfDeltaVariant = (dfVariant1 - dfVariant) /
349 1 : (double)nDeltaY;
350 :
351 33 : while ( nDeltaY-- >= 0 )
352 : {
353 31 : if ( 0 <= iX && iX < nRasterXSize
354 : && 0 <= iY && iY < nRasterYSize )
355 31 : pfnPointFunc( pCBData, iY, iX, dfVariant );
356 :
357 31 : dfVariant += dfDeltaVariant;
358 31 : iY += nYStep;
359 31 : if ( nError > 0 )
360 : {
361 15 : iX += nXStep;
362 15 : nError += nYError;
363 : }
364 : else
365 16 : nError += nXError;
366 : }
367 : }
368 : }
369 : }
370 : }
371 :
372 : /************************************************************************/
373 : /* GDALdllImageLineAllTouched() */
374 : /* */
375 : /* This alternate line drawing algorithm attempts to ensure */
376 : /* that every pixel touched at all by the line will get set. */
377 : /* @param padfVariant should contain the values that are to be */
378 : /* added to the burn value. The values along the line between the */
379 : /* points will be linearly interpolated. These values are used */
380 : /* only if pCBData->eBurnValueSource is set to something other */
381 : /* than GBV_UserBurnValue. If NULL is passed, a monotonous line */
382 : /* will be drawn with the burn value. */
383 : /************************************************************************/
384 :
385 : void
386 7 : GDALdllImageLineAllTouched(int nRasterXSize, int nRasterYSize,
387 : int nPartCount, int *panPartSize,
388 : double *padfX, double *padfY, double *padfVariant,
389 : llPointFunc pfnPointFunc, void *pCBData )
390 :
391 : {
392 : int i, n;
393 :
394 7 : if ( !nPartCount )
395 0 : return;
396 :
397 14 : for ( i = 0, n = 0; i < nPartCount; n += panPartSize[i++] )
398 : {
399 : int j;
400 :
401 33 : for ( j = 1; j < panPartSize[i]; j++ )
402 : {
403 26 : double dfX = padfX[n + j - 1];
404 26 : double dfY = padfY[n + j - 1];
405 :
406 26 : double dfXEnd = padfX[n + j];
407 26 : double dfYEnd = padfY[n + j];
408 :
409 26 : double dfVariant = 0, dfVariantEnd = 0;
410 26 : if( padfVariant != NULL &&
411 : ((GDALRasterizeInfo *)pCBData)->eBurnValueSource !=
412 : GBV_UserBurnValue )
413 : {
414 0 : dfVariant = padfVariant[n + j - 1];
415 0 : dfVariantEnd = padfVariant[n + j];
416 : }
417 :
418 : // Skip segments that are off the target region.
419 26 : if( (dfY < 0 && dfYEnd < 0)
420 : || (dfY > nRasterYSize && dfYEnd > nRasterYSize)
421 : || (dfX < 0 && dfXEnd < 0)
422 : || (dfX > nRasterXSize && dfXEnd > nRasterXSize) )
423 0 : continue;
424 :
425 : // Swap if needed so we can proceed from left2right (X increasing)
426 26 : if( dfX > dfXEnd )
427 : {
428 9 : llSwapDouble( &dfX, &dfXEnd );
429 9 : llSwapDouble( &dfY, &dfYEnd );
430 9 : llSwapDouble( &dfVariant, &dfVariantEnd );
431 : }
432 :
433 : // Special case for vertical lines.
434 26 : if( floor(dfX) == floor(dfXEnd) )
435 : {
436 12 : if( dfYEnd < dfY )
437 : {
438 7 : llSwapDouble( &dfY, &dfYEnd );
439 7 : llSwapDouble( &dfVariant, &dfVariantEnd );
440 : }
441 :
442 12 : int iX = (int) floor(dfX);
443 12 : int iY = (int) floor(dfY);
444 12 : int iYEnd = (int) floor(dfYEnd);
445 :
446 12 : if( iX >= nRasterXSize )
447 0 : continue;
448 :
449 12 : double dfDeltaVariant = 0;
450 12 : if(( dfYEnd - dfY ) > 0)
451 : dfDeltaVariant = ( dfVariantEnd - dfVariant )
452 12 : / ( dfYEnd - dfY );//per unit change in iY
453 :
454 : // Clip to the borders of the target region
455 12 : if( iY < 0 )
456 0 : iY = 0;
457 12 : if( iYEnd >= nRasterYSize )
458 0 : iYEnd = nRasterYSize - 1;
459 12 : dfVariant += dfDeltaVariant * ( (double)iY - dfY );
460 :
461 12 : if( padfVariant == NULL )
462 112 : for( ; iY <= iYEnd; iY++ )
463 100 : pfnPointFunc( pCBData, iY, iX, 0.0 );
464 : else
465 0 : for( ; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant )
466 0 : pfnPointFunc( pCBData, iY, iX, dfVariant );
467 :
468 12 : continue; // next segment
469 : }
470 :
471 : double dfDeltaVariant = ( dfVariantEnd - dfVariant )
472 14 : / ( dfXEnd - dfX );//per unit change in iX
473 :
474 : // special case for horizontal lines
475 14 : if( floor(dfY) == floor(dfYEnd) )
476 : {
477 12 : if( dfXEnd < dfX )
478 : {
479 0 : llSwapDouble( &dfX, &dfXEnd );
480 0 : llSwapDouble( &dfVariant, &dfVariantEnd );
481 : }
482 :
483 12 : int iX = (int) floor(dfX);
484 12 : int iY = (int) floor(dfY);
485 12 : int iXEnd = (int) floor(dfXEnd);
486 :
487 12 : if( iY >= nRasterYSize )
488 0 : continue;
489 :
490 : // Clip to the borders of the target region
491 12 : if( iX < 0 )
492 0 : iX = 0;
493 12 : if( iXEnd >= nRasterXSize )
494 0 : iXEnd = nRasterXSize - 1;
495 12 : dfVariant += dfDeltaVariant * ( (double)iX - dfX );
496 :
497 12 : if( padfVariant == NULL )
498 136 : for( ; iX <= iXEnd; iX++ )
499 124 : pfnPointFunc( pCBData, iY, iX, 0.0 );
500 : else
501 0 : for( ; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant )
502 0 : pfnPointFunc( pCBData, iY, iX, dfVariant );
503 :
504 12 : continue; // next segment
505 : }
506 :
507 : /* -------------------------------------------------------------------- */
508 : /* General case - left to right sloped. */
509 : /* -------------------------------------------------------------------- */
510 2 : double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
511 :
512 : // clip segment in X
513 2 : if( dfXEnd > nRasterXSize )
514 : {
515 0 : dfYEnd -= ( dfXEnd - (double)nRasterXSize ) * dfSlope;
516 0 : dfXEnd = nRasterXSize;
517 : }
518 2 : if( dfX < 0 )
519 : {
520 0 : dfY += (0 - dfX) * dfSlope;
521 0 : dfVariant += dfDeltaVariant * (0.0 - dfX);
522 0 : dfX = 0.0;
523 : }
524 :
525 : // clip segment in Y
526 : double dfDiffX;
527 2 : if( dfYEnd > dfY )
528 : {
529 0 : if( dfY < 0 )
530 : {
531 0 : dfX += (dfDiffX = (0 - dfY) / dfSlope);
532 0 : dfVariant += dfDeltaVariant * dfDiffX;
533 0 : dfY = 0.0;
534 : }
535 0 : if( dfYEnd >= nRasterYSize )
536 : {
537 0 : dfXEnd += ( dfYEnd - (double)nRasterYSize ) / dfSlope;
538 0 : dfYEnd = nRasterXSize;
539 : }
540 : }
541 : else
542 : {
543 2 : if( dfY >= nRasterYSize )
544 : {
545 0 : dfX += (dfDiffX = ((double)nRasterYSize - dfY) / dfSlope);
546 0 : dfVariant += dfDeltaVariant * dfDiffX;
547 0 : dfY = nRasterYSize;
548 : }
549 2 : if( dfYEnd < 0 )
550 : {
551 0 : dfXEnd -= ( dfYEnd - 0 ) / dfSlope;
552 0 : dfYEnd = 0;
553 : }
554 : }
555 :
556 : // step from pixel to pixel.
557 14 : while( dfX < dfXEnd )
558 : {
559 10 : int iX = (int) floor(dfX);
560 10 : int iY = (int) floor(dfY);
561 :
562 : // burn in the current point.
563 : // We should be able to drop the Y check because we cliped in Y,
564 : // but there may be some error with all the small steps.
565 10 : if( iY >= 0 && iY < nRasterYSize )
566 10 : pfnPointFunc( pCBData, iY, iX, dfVariant );
567 :
568 10 : double dfStepX = floor(dfX+1.0) - dfX;
569 10 : double dfStepY = dfStepX * dfSlope;
570 :
571 : // step to right pixel without changing scanline?
572 10 : if( (int) floor(dfY + dfStepY) == iY )
573 : {
574 6 : dfX += dfStepX;
575 6 : dfY += dfStepY;
576 6 : dfVariant += dfDeltaVariant * dfStepX;
577 : }
578 4 : else if( dfSlope < 0 )
579 : {
580 4 : dfStepY = iY - dfY;
581 4 : if( dfStepY > -0.000000001 )
582 4 : dfStepY = -0.000000001;
583 :
584 4 : dfStepX = dfStepY / dfSlope;
585 4 : dfX += dfStepX;
586 4 : dfY += dfStepY;
587 4 : dfVariant += dfDeltaVariant * dfStepX;
588 : }
589 : else
590 : {
591 0 : dfStepY = (iY + 1) - dfY;
592 0 : if( dfStepY < 0.000000001 )
593 0 : dfStepY = 0.000000001;
594 :
595 0 : dfStepX = dfStepY / dfSlope;
596 0 : dfX += dfStepX;
597 0 : dfY += dfStepY;
598 0 : dfVariant += dfDeltaVariant * dfStepX;
599 : }
600 : } // next step along sement.
601 :
602 : } // next segment
603 : } // next part
604 : }
605 :
|
