PLplot 5.9.6
|
00001 /* $Id$ 00002 * 00003 * Polygon pattern fill. 00004 * 00005 * Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 Alan W. Irwin 00006 * Copyright (C) 2005, 2006, 2007, 2008, 2009 Arjen Markus 00007 * 00008 * This file is part of PLplot. 00009 * 00010 * PLplot is free software; you can redistribute it and/or modify 00011 * it under the terms of the GNU General Library Public License as published 00012 * by the Free Software Foundation; either version 2 of the License, or 00013 * (at your option) any later version. 00014 * 00015 * PLplot is distributed in the hope that it will be useful, 00016 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 * GNU Library General Public License for more details. 00019 * 00020 * You should have received a copy of the GNU Library General Public License 00021 * along with PLplot; if not, write to the Free Software 00022 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00023 * 00024 */ 00025 00026 #include "plplotP.h" 00027 00028 #define INSIDE( ix, iy ) ( BETW( ix, xmin, xmax ) && BETW( iy, ymin, ymax ) ) 00029 00030 #define DTOR ( PI / 180. ) 00031 #define BINC 50 00032 /* Near-border comparison criterion (NBCC). */ 00033 #define PL_NBCC 2 00034 /* Variant of BETW that returns true if between or within PL_NBCC of it. */ 00035 #define BETW_NBCC( ix, ia, ib ) ( ( ( ix ) <= ( ia + PL_NBCC ) && ( ix ) >= ( ib - PL_NBCC ) ) || ( ( ix ) >= ( ia - PL_NBCC ) && ( ix ) <= ( ib + PL_NBCC ) ) ) 00036 00037 /* Status codes ORed together in the return from notcrossed. 00038 * PL_NOT_CROSSED is set whenever the two line segments definitely 00039 * (i.e., intersection not near the ends or completely apart) 00040 * do not cross each other. 00041 * 00042 * PL_NEAR_A1 is set if the intersect is near (+/- PL_NBCC) the first 00043 * end of line segment A. 00044 * 00045 * PL_NEAR_A2 is set if the intersect is near (+/- PL_NBCC) the second 00046 * end of line segment A. 00047 * 00048 * PL_NEAR_B1 is set if the intersect is near (+/- PL_NBCC) the first 00049 * end of line segment B. 00050 * 00051 * PL_NEAR_B2 is set if the intersect is near (+/- PL_NBCC) the second 00052 * end of line segment B. 00053 * 00054 * PL_NEAR_PARALLEL is set if the two line segments are nearly parallel 00055 * with each other, i.e., a change in one of the coordinates of up to 00056 * (+/- PL_NBCC) would render them exactly parallel. 00057 * 00058 * PL_PARALLEL is set if the two line segments are exactly parallel 00059 * with each other. 00060 */ 00061 enum PL_CrossedStatus 00062 { 00063 PL_NOT_CROSSED = 0x1, 00064 PL_NEAR_A1 = 0x2, 00065 PL_NEAR_A2 = 0x4, 00066 PL_NEAR_B1 = 0x8, 00067 PL_NEAR_B2 = 0x10, 00068 PL_NEAR_PARALLEL = 0x20, 00069 PL_PARALLEL = 0x40 00070 }; 00071 00072 struct point 00073 { 00074 PLINT x, y; 00075 }; 00076 static PLINT bufferleng, buffersize, *buffer; 00077 00078 /* Static function prototypes */ 00079 00080 static int 00081 compar( const void *, const void * ); 00082 00083 static void 00084 addcoord( PLINT, PLINT ); 00085 00086 static void 00087 tran( PLINT *, PLINT *, PLFLT, PLFLT ); 00088 00089 static void 00090 buildlist( PLINT, PLINT, PLINT, PLINT, PLINT, PLINT, PLINT ); 00091 00092 static int 00093 notpointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ); 00094 00095 static int 00096 circulation( PLINT *x, PLINT *y, PLINT npts ); 00097 00098 static void 00099 fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, 00100 PLINT fill_status, 00101 void ( *fill )( short *, short *, PLINT ), 00102 const PLINT *x1, const PLINT *y1, 00103 PLINT i1start, PLINT n1, 00104 const PLINT *x2, const PLINT *y2, 00105 const PLINT *if2, PLINT n2 ); 00106 00107 static int 00108 notcrossed( PLINT *xintersect, PLINT *yintersect, 00109 PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, 00110 PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); 00111 00112 static int 00113 positive_orientation( PLINT n, const PLINT *x, const PLINT *y ); 00114 00115 static int 00116 number_crossings( PLINT *xcross, PLINT *ycross, PLINT *i2cross, PLINT ncross, 00117 PLINT i1, PLINT n1, const PLINT *x1, const PLINT *y1, 00118 PLINT n2, const PLINT *x2, const PLINT *y2 ); 00119 00120 00121 /*----------------------------------------------------------------------*\ 00122 * void plfill() 00123 * 00124 * Pattern fills the polygon bounded by the input points. 00125 * If hardware fill is used, a maximum of PL_MAXPOLY-1 vertices is allowed. 00126 * The final point is explicitly added if it doesn't match up to the first, 00127 * to prevent clipping problems. 00128 \*----------------------------------------------------------------------*/ 00129 00130 void 00131 c_plfill( PLINT n, PLFLT *x, PLFLT *y ) 00132 { 00133 PLINT xpoly[PL_MAXPOLY], ypoly[PL_MAXPOLY]; 00134 PLINT i; 00135 PLFLT xt, yt; 00136 00137 if ( plsc->level < 3 ) 00138 { 00139 plabort( "plfill: Please set up window first" ); 00140 return; 00141 } 00142 if ( n < 3 ) 00143 { 00144 plabort( "plfill: Not enough points in object" ); 00145 return; 00146 } 00147 if ( n > PL_MAXPOLY - 1 ) 00148 { 00149 plwarn( "plfill: too many points in polygon" ); 00150 n = PL_MAXPOLY; 00151 } 00152 for ( i = 0; i < n; i++ ) 00153 { 00154 TRANSFORM( x[i], y[i], &xt, &yt ); 00155 xpoly[i] = plP_wcpcx( xt ); 00156 ypoly[i] = plP_wcpcy( yt ); 00157 } 00158 00159 if ( x[0] != x[n - 1] || y[0] != y[n - 1] ) 00160 { 00161 if ( n < PL_MAXPOLY ) 00162 n++; 00163 TRANSFORM( x[0], y[0], &xt, &yt ); 00164 xpoly[n - 1] = plP_wcpcx( xt ); 00165 ypoly[n - 1] = plP_wcpcy( yt ); 00166 } 00167 00168 plP_plfclp( xpoly, ypoly, n, plsc->clpxmi, plsc->clpxma, 00169 plsc->clpymi, plsc->clpyma, plP_fill ); 00170 } 00171 00172 /*----------------------------------------------------------------------*\ 00173 * void plfill3() 00174 * 00175 * Pattern fills the polygon in 3d bounded by the input points. 00176 * If hardware fill is used, a maximum of PL_MAXPOLY-1 vertices is allowed. 00177 * The final point is explicitly added if it doesn't match up to the first, 00178 * to prevent clipping problems. 00179 \*----------------------------------------------------------------------*/ 00180 00181 void 00182 c_plfill3( PLINT n, PLFLT *x, PLFLT *y, PLFLT *z ) 00183 { 00184 PLFLT tx[PL_MAXPOLY], ty[PL_MAXPOLY], tz[PL_MAXPOLY]; 00185 PLFLT *V[3]; 00186 PLINT xpoly[PL_MAXPOLY], ypoly[PL_MAXPOLY]; 00187 PLINT i; 00188 PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale; 00189 00190 if ( plsc->level < 3 ) 00191 { 00192 plabort( "plfill3: Please set up window first" ); 00193 return; 00194 } 00195 if ( n < 3 ) 00196 { 00197 plabort( "plfill3: Not enough points in object" ); 00198 return; 00199 } 00200 if ( n > PL_MAXPOLY - 1 ) 00201 { 00202 plwarn( "plfill3: too many points in polygon" ); 00203 n = PL_MAXPOLY; 00204 } 00205 00206 plP_gdom( &xmin, &xmax, &ymin, &ymax ); 00207 plP_grange( &zscale, &zmin, &zmax ); 00208 00209 /* copy the vertices so we can clip without corrupting the input */ 00210 for ( i = 0; i < n; i++ ) 00211 { 00212 tx[i] = x[i]; ty[i] = y[i]; tz[i] = z[i]; 00213 } 00214 if ( tx[0] != tx[n - 1] || ty[0] != ty[n - 1] || tz[0] != tz[n - 1] ) 00215 { 00216 if ( n < PL_MAXPOLY ) 00217 n++; 00218 tx[n - 1] = tx[0]; ty[n - 1] = ty[0]; tz[n - 1] = tz[0]; 00219 } 00220 V[0] = tx; V[1] = ty; V[2] = tz; 00221 n = plP_clip_poly( n, V, 0, 1, -xmin ); 00222 n = plP_clip_poly( n, V, 0, -1, xmax ); 00223 n = plP_clip_poly( n, V, 1, 1, -ymin ); 00224 n = plP_clip_poly( n, V, 1, -1, ymax ); 00225 n = plP_clip_poly( n, V, 2, 1, -zmin ); 00226 n = plP_clip_poly( n, V, 2, -1, zmax ); 00227 for ( i = 0; i < n; i++ ) 00228 { 00229 xpoly[i] = plP_wcpcx( plP_w3wcx( tx[i], ty[i], tz[i] ) ); 00230 ypoly[i] = plP_wcpcy( plP_w3wcy( tx[i], ty[i], tz[i] ) ); 00231 } 00232 00233 /* AWI: in the past we have used 00234 * plP_fill(xpoly, ypoly, n); 00235 * here, but our educated guess is this fill should be done via the clipping 00236 * interface instead as below. 00237 * No example tests this code so one of our users will end up inadvertently 00238 * testing this for us. 00239 * 00240 * jc: I have checked, and both versions does give the same result, i.e., clipping 00241 * to the window boundaries. The reason is that the above plP_clip_poly() does 00242 * the clipping. To check this, is enough to diminish the x/y/z min/max arguments in 00243 * plw3d() in x08c. But let's keep it, although 10% slower... 00244 */ 00245 plP_plfclp( xpoly, ypoly, n, plsc->clpxmi, plsc->clpxma, 00246 plsc->clpymi, plsc->clpyma, plP_fill ); 00247 } 00248 00249 /*----------------------------------------------------------------------*\ 00250 * void plfill_soft() 00251 * 00252 * Pattern fills in software the polygon bounded by the input points. 00253 \*----------------------------------------------------------------------*/ 00254 00255 void 00256 plfill_soft( short *x, short *y, PLINT n ) 00257 { 00258 PLINT i, j; 00259 PLINT xp1, yp1, xp2, yp2, xp3, yp3; 00260 PLINT k, dinc; 00261 PLFLT ci, si; 00262 double temp; 00263 00264 buffersize = 2 * BINC; 00265 buffer = (PLINT *) malloc( (size_t) buffersize * sizeof ( PLINT ) ); 00266 if ( !buffer ) 00267 { 00268 plabort( "plfill: Out of memory" ); 00269 return; 00270 } 00271 00272 /* Loop over sets of lines in pattern */ 00273 00274 for ( k = 0; k < plsc->nps; k++ ) 00275 { 00276 bufferleng = 0; 00277 00278 temp = DTOR * plsc->inclin[k] * 0.1; 00279 si = sin( temp ) * plsc->ypmm; 00280 ci = cos( temp ) * plsc->xpmm; 00281 00282 /* normalize: 1 = si*si + ci*ci */ 00283 00284 temp = sqrt( (double) ( si * si + ci * ci ) ); 00285 si /= temp; 00286 ci /= temp; 00287 00288 dinc = (PLINT) ( plsc->delta[k] * SSQR( plsc->ypmm * ABS( ci ), 00289 plsc->xpmm * ABS( si ) ) / 1000. ); 00290 00291 if ( dinc < 0 ) 00292 dinc = -dinc; 00293 if ( dinc == 0 ) 00294 dinc = 1; 00295 00296 xp1 = x[n - 2]; 00297 yp1 = y[n - 2]; 00298 tran( &xp1, &yp1, (PLFLT) ci, (PLFLT) si ); 00299 00300 xp2 = x[n - 1]; 00301 yp2 = y[n - 1]; 00302 tran( &xp2, &yp2, (PLFLT) ci, (PLFLT) si ); 00303 00304 /* Loop over points in polygon */ 00305 00306 for ( i = 0; i < n; i++ ) 00307 { 00308 xp3 = x[i]; 00309 yp3 = y[i]; 00310 tran( &xp3, &yp3, (PLFLT) ci, (PLFLT) si ); 00311 buildlist( xp1, yp1, xp2, yp2, xp3, yp3, dinc ); 00312 xp1 = xp2; 00313 yp1 = yp2; 00314 xp2 = xp3; 00315 yp2 = yp3; 00316 } 00317 00318 /* Sort list by y then x */ 00319 00320 qsort( (void *) buffer, (size_t) bufferleng / 2, 00321 (size_t) sizeof ( struct point ), compar ); 00322 00323 /* OK, now do the hatching */ 00324 00325 i = 0; 00326 00327 while ( i < bufferleng ) 00328 { 00329 xp1 = buffer[i]; 00330 yp1 = buffer[i + 1]; 00331 i += 2; 00332 xp2 = xp1; 00333 yp2 = yp1; 00334 tran( &xp1, &yp1, (PLFLT) ci, (PLFLT) ( -si ) ); 00335 plP_movphy( xp1, yp1 ); 00336 xp1 = buffer[i]; 00337 yp1 = buffer[i + 1]; 00338 i += 2; 00339 if ( yp2 != yp1 ) 00340 { 00341 fprintf( stderr, "plfill: oh oh we are lost\n" ); 00342 for ( j = 0; j < bufferleng; j += 2 ) 00343 { 00344 fprintf( stderr, "plfill: %d %d\n", 00345 (int) buffer[j], (int) buffer[j + 1] ); 00346 } 00347 continue; /* Uh oh we're lost */ 00348 } 00349 tran( &xp1, &yp1, (PLFLT) ci, (PLFLT) ( -si ) ); 00350 plP_draphy( xp1, yp1 ); 00351 } 00352 } 00353 free( (void *) buffer ); 00354 } 00355 00356 /*----------------------------------------------------------------------*\ 00357 * Utility functions 00358 \*----------------------------------------------------------------------*/ 00359 00360 void 00361 tran( PLINT *a, PLINT *b, PLFLT c, PLFLT d ) 00362 { 00363 PLINT ta, tb; 00364 00365 ta = *a; 00366 tb = *b; 00367 00368 *a = (PLINT) floor( (double) ( ta * c + tb * d + 0.5 ) ); 00369 *b = (PLINT) floor( (double) ( tb * c - ta * d + 0.5 ) ); 00370 } 00371 00372 void 00373 buildlist( PLINT xp1, PLINT yp1, PLINT xp2, PLINT yp2, PLINT xp3, PLINT yp3, 00374 PLINT dinc ) 00375 { 00376 PLINT min_y, max_y; 00377 PLINT dx, dy, cstep, nstep, ploty, plotx; 00378 00379 dx = xp2 - xp1; 00380 dy = yp2 - yp1; 00381 00382 if ( dy == 0 ) 00383 { 00384 if ( yp2 > yp3 && ( ( yp2 % dinc ) == 0 ) ) 00385 addcoord( xp2, yp2 ); 00386 return; 00387 } 00388 00389 if ( dy > 0 ) 00390 { 00391 cstep = 1; 00392 min_y = yp1; 00393 max_y = yp2; 00394 } 00395 else 00396 { 00397 cstep = -1; 00398 min_y = yp2; 00399 max_y = yp1; 00400 } 00401 00402 nstep = ( yp3 > yp2 ? 1 : -1 ); 00403 if ( yp3 == yp2 ) 00404 nstep = 0; 00405 00406 /* Build coordinate list */ 00407 00408 ploty = ( min_y / dinc ) * dinc; 00409 if ( ploty < min_y ) 00410 ploty += dinc; 00411 00412 for (; ploty <= max_y; ploty += dinc ) 00413 { 00414 if ( ploty == yp1 ) 00415 continue; 00416 if ( ploty == yp2 ) 00417 { 00418 if ( cstep == -nstep ) 00419 continue; 00420 if ( yp2 == yp3 && yp1 > yp2 ) 00421 continue; 00422 } 00423 plotx = xp1 + (PLINT) floor( ( (double) ( ploty - yp1 ) * dx ) / dy + 0.5 ); 00424 addcoord( plotx, ploty ); 00425 } 00426 } 00427 00428 void 00429 addcoord( PLINT xp1, PLINT yp1 ) 00430 { 00431 PLINT *temp; 00432 00433 if ( bufferleng + 2 > buffersize ) 00434 { 00435 buffersize += 2 * BINC; 00436 temp = (PLINT *) realloc( (void *) buffer, 00437 (size_t) buffersize * sizeof ( PLINT ) ); 00438 if ( !temp ) 00439 { 00440 free( (void *) buffer ); 00441 plexit( "plfill: Out of memory!" ); 00442 } 00443 buffer = temp; 00444 } 00445 00446 buffer[bufferleng++] = xp1; 00447 buffer[bufferleng++] = yp1; 00448 } 00449 00450 int 00451 compar( const void *pnum1, const void *pnum2 ) 00452 { 00453 const struct point *pnt1, *pnt2; 00454 00455 pnt1 = (const struct point *) pnum1; 00456 pnt2 = (const struct point *) pnum2; 00457 00458 if ( pnt1->y < pnt2->y ) 00459 return -1; 00460 else if ( pnt1->y > pnt2->y ) 00461 return 1; 00462 00463 /* Only reach here if y coords are equal, so sort by x */ 00464 00465 if ( pnt1->x < pnt2->x ) 00466 return -1; 00467 else if ( pnt1->x > pnt2->x ) 00468 return 1; 00469 00470 return 0; 00471 } 00472 00473 /*----------------------------------------------------------------------*\ 00474 * void plP_plfclp() 00475 * 00476 * Fills a polygon within the clip limits. 00477 \*----------------------------------------------------------------------*/ 00478 00479 void 00480 plP_plfclp( PLINT *x, PLINT *y, PLINT npts, 00481 PLINT xmin, PLINT xmax, PLINT ymin, PLINT ymax, 00482 void ( *draw )( short *, short *, PLINT ) ) 00483 { 00484 #ifdef USE_FILL_INTERSECTION_POLYGON 00485 PLINT *x10, *y10, *x1, *y1, *if1, i1start = 0, i, im1, n1, n1m1, 00486 ifnotpointinpolygon; 00487 PLINT x2[4] = { xmin, xmax, xmax, xmin }; 00488 PLINT y2[4] = { ymin, ymin, ymax, ymax }; 00489 PLINT if2[4] = { 0, 0, 0, 0 }; 00490 PLINT n2 = 4; 00491 00492 /* Must have at least 3 points and draw() specified */ 00493 if ( npts < 3 || !draw ) 00494 return; 00495 00496 if ( ( x10 = (PLINT *) malloc( npts * sizeof ( PLINT ) ) ) == NULL ) 00497 { 00498 plexit( "plP_plfclp: Insufficient memory" ); 00499 } 00500 if ( ( y10 = (PLINT *) malloc( npts * sizeof ( PLINT ) ) ) == NULL ) 00501 { 00502 plexit( "plP_plfclp: Insufficient memory" ); 00503 } 00504 /* Polygon 2 obviously has no dups nor two consective segments that 00505 * are parallel, but get rid of those type of segments in polygon 1 00506 * if they exist. */ 00507 00508 im1 = npts - 1; 00509 n1 = 0; 00510 for ( i = 0; i < npts; i++ ) 00511 { 00512 if ( !( x[i] == x[im1] && y[i] == y[im1] ) ) 00513 { 00514 x10[n1] = x[i]; 00515 y10[n1++] = y[i]; 00516 } 00517 im1 = i; 00518 } 00519 00520 /* Must have at least three points that satisfy the above criteria. */ 00521 if ( n1 < 3 ) 00522 { 00523 free( x10 ); 00524 free( y10 ); 00525 return; 00526 } 00527 00528 /* Polygon 2 obviously has a positive orientation (i.e., as you 00529 * ascend in index along the boundary, the points just adjacent to 00530 * the boundary and on the left are interior points for the 00531 * polygon), but enforce this condition demanded by 00532 * fill_intersection_polygon for polygon 1 as well. */ 00533 if ( positive_orientation( n1, x10, y10 ) ) 00534 { 00535 x1 = x10; 00536 y1 = y10; 00537 } 00538 else 00539 { 00540 if ( ( x1 = (PLINT *) malloc( n1 * sizeof ( PLINT ) ) ) == NULL ) 00541 { 00542 plexit( "plP_plfclp: Insufficient memory" ); 00543 } 00544 if ( ( y1 = (PLINT *) malloc( n1 * sizeof ( PLINT ) ) ) == NULL ) 00545 { 00546 plexit( "plP_plfclp: Insufficient memory" ); 00547 } 00548 n1m1 = n1 - 1; 00549 for ( i = 0; i < n1; i++ ) 00550 { 00551 x1[n1m1 - i] = x10[i]; 00552 y1[n1m1 - i] = y10[i]; 00553 } 00554 free( x10 ); 00555 free( y10 ); 00556 } 00557 00558 /* Insure that the first vertex of polygon 1 (starting with n1 - 00559 * 1) that is not on the border of polygon 2 is definitely outside 00560 * polygon 2. */ 00561 im1 = n1 - 1; 00562 for ( i = 0; i < n1; i++ ) 00563 { 00564 if ( ( ifnotpointinpolygon = 00565 notpointinpolygon( n2, x2, y2, x1[im1], y1[im1] ) ) != 1 ) 00566 break; 00567 im1 = i; 00568 } 00569 00570 if ( ifnotpointinpolygon ) 00571 fill_intersection_polygon( 0, 0, 0, draw, x1, y1, i1start, n1, x2, y2, if2, n2 ); 00572 else 00573 { 00574 if ( ( if1 = (PLINT *) calloc( n1, sizeof ( PLINT ) ) ) == NULL ) 00575 { 00576 plexit( "plP_plfclp: Insufficient memory" ); 00577 } 00578 fill_intersection_polygon( 0, 0, 0, draw, x2, y2, i1start, n2, x1, y1, if1, n1 ); 00579 free( if1 ); 00580 } 00581 free( x1 ); 00582 free( y1 ); 00583 return; 00584 } 00585 #else /* USE_FILL_INTERSECTION_POLYGON */ 00586 00587 PLINT i, x1, x2, y1, y2; 00588 int iclp = 0, iout = 2; 00589 short _xclp[2 * PL_MAXPOLY + 2], _yclp[2 * PL_MAXPOLY + 2]; 00590 short *xclp, *yclp; 00591 int drawable; 00592 int crossed_xmin1 = 0, crossed_xmax1 = 0; 00593 int crossed_ymin1 = 0, crossed_ymax1 = 0; 00594 int crossed_xmin2 = 0, crossed_xmax2 = 0; 00595 int crossed_ymin2 = 0, crossed_ymax2 = 0; 00596 int crossed_up = 0, crossed_down = 0; 00597 int crossed_left = 0, crossed_right = 0; 00598 int inside_lb; 00599 int inside_lu; 00600 int inside_rb; 00601 int inside_ru; 00602 00603 /* Must have at least 3 points and draw() specified */ 00604 if ( npts < 3 || !draw ) 00605 return; 00606 00607 if ( npts < PL_MAXPOLY ) 00608 { 00609 xclp = _xclp; 00610 yclp = _yclp; 00611 } 00612 else 00613 { 00614 if ( ( ( xclp = (short *) malloc( ( 2 * npts + 2 ) * sizeof ( short ) ) ) == NULL ) || 00615 ( ( yclp = (short *) malloc( ( 2 * npts + 2 ) * sizeof ( short ) ) ) == NULL ) ) 00616 { 00617 plexit( "plP_plfclp: Insufficient memory" ); 00618 } 00619 } 00620 inside_lb = !notpointinpolygon( npts, x, y, xmin, ymin ); 00621 inside_lu = !notpointinpolygon( npts, x, y, xmin, ymax ); 00622 inside_rb = !notpointinpolygon( npts, x, y, xmax, ymin ); 00623 inside_ru = !notpointinpolygon( npts, x, y, xmax, ymax ); 00624 00625 for ( i = 0; i < npts - 1; i++ ) 00626 { 00627 x1 = x[i]; x2 = x[i + 1]; 00628 y1 = y[i]; y2 = y[i + 1]; 00629 00630 drawable = ( INSIDE( x1, y1 ) && INSIDE( x2, y2 ) ); 00631 if ( !drawable ) 00632 drawable = !plP_clipline( &x1, &y1, &x2, &y2, 00633 xmin, xmax, ymin, ymax ); 00634 00635 if ( drawable ) 00636 { 00637 /* Boundary crossing condition -- coming in. */ 00638 crossed_xmin2 = ( x1 == xmin ); crossed_xmax2 = ( x1 == xmax ); 00639 crossed_ymin2 = ( y1 == ymin ); crossed_ymax2 = ( y1 == ymax ); 00640 00641 crossed_left = ( crossed_left || crossed_xmin2 ); 00642 crossed_right = ( crossed_right || crossed_xmax2 ); 00643 crossed_down = ( crossed_down || crossed_ymin2 ); 00644 crossed_up = ( crossed_up || crossed_ymax2 ); 00645 iout = iclp + 2; 00646 /* If the first segment, just add it. */ 00647 00648 if ( iclp == 0 ) 00649 { 00650 xclp[iclp] = x1; yclp[iclp] = y1; iclp++; 00651 xclp[iclp] = x2; yclp[iclp] = y2; iclp++; 00652 } 00653 00654 /* Not first point. If first point of this segment matches up to the 00655 * previous point, just add it. */ 00656 00657 else if ( x1 == xclp[iclp - 1] && y1 == yclp[iclp - 1] ) 00658 { 00659 xclp[iclp] = x2; yclp[iclp] = y2; iclp++; 00660 } 00661 00662 /* Otherwise, we need to add both points, to connect the points in the 00663 * polygon along the clip boundary. If we encircled a corner, we have 00664 * to add that first. 00665 */ 00666 00667 else 00668 { 00669 /* Treat the case where we encircled two corners: 00670 * Construct a polygon out of the subset of vertices 00671 * Note that the direction is important too when adding 00672 * the extra points */ 00673 xclp[iclp + 1] = x2; yclp[iclp + 1] = y2; 00674 xclp[iclp + 2] = x1; yclp[iclp + 2] = y1; 00675 iout = iout - iclp + 1; 00676 /* Upper two */ 00677 if ( ( ( crossed_xmin1 && crossed_xmax2 ) || 00678 ( crossed_xmin2 && crossed_xmax1 ) ) && 00679 inside_lu ) 00680 { 00681 if ( crossed_xmin1 ) 00682 { 00683 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00684 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00685 } 00686 else 00687 { 00688 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00689 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00690 } 00691 } 00692 /* Lower two */ 00693 else if ( ( ( crossed_xmin1 && crossed_xmax2 ) || 00694 ( crossed_xmin2 && crossed_xmax1 ) ) && 00695 inside_lb ) 00696 { 00697 if ( crossed_xmin1 ) 00698 { 00699 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00700 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00701 } 00702 else 00703 { 00704 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00705 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00706 } 00707 } 00708 /* Left two */ 00709 else if ( ( ( crossed_ymin1 && crossed_ymax2 ) || 00710 ( crossed_ymin2 && crossed_ymax1 ) ) && 00711 inside_lb ) 00712 { 00713 if ( crossed_ymin1 ) 00714 { 00715 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00716 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00717 } 00718 else 00719 { 00720 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00721 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00722 } 00723 } 00724 /* Right two */ 00725 else if ( ( ( crossed_ymin1 && crossed_ymax2 ) || 00726 ( crossed_ymin2 && crossed_ymax1 ) ) && 00727 inside_rb ) 00728 { 00729 if ( crossed_ymin1 ) 00730 { 00731 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00732 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00733 } 00734 else 00735 { 00736 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00737 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00738 } 00739 } 00740 /* Now the case where we encircled one corner */ 00741 /* Lower left */ 00742 else if ( ( crossed_xmin1 && crossed_ymin2 ) || 00743 ( crossed_ymin1 && crossed_xmin2 ) ) 00744 { 00745 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00746 } 00747 /* Lower right */ 00748 else if ( ( crossed_xmax1 && crossed_ymin2 ) || 00749 ( crossed_ymin1 && crossed_xmax2 ) ) 00750 { 00751 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00752 } 00753 /* Upper left */ 00754 else if ( ( crossed_xmin1 && crossed_ymax2 ) || 00755 ( crossed_ymax1 && crossed_xmin2 ) ) 00756 { 00757 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00758 } 00759 /* Upper right */ 00760 else if ( ( crossed_xmax1 && crossed_ymax2 ) || 00761 ( crossed_ymax1 && crossed_xmax2 ) ) 00762 { 00763 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00764 } 00765 00766 /* Now add current segment. */ 00767 xclp[iclp] = x1; yclp[iclp] = y1; iclp++; 00768 xclp[iclp] = x2; yclp[iclp] = y2; iclp++; 00769 } 00770 00771 /* Boundary crossing condition -- going out. */ 00772 crossed_xmin1 = ( x2 == xmin ); crossed_xmax1 = ( x2 == xmax ); 00773 crossed_ymin1 = ( y2 == ymin ); crossed_ymax1 = ( y2 == ymax ); 00774 } 00775 } 00776 00777 /* Limit case - all vertices are outside of bounding box. So just fill entire 00778 * box, *if* the bounding box is completely encircled. 00779 */ 00780 if ( iclp == 0 ) 00781 { 00782 if ( inside_lb ) 00783 { 00784 xclp[0] = xmin; yclp[0] = ymin; 00785 xclp[1] = xmax; yclp[1] = ymin; 00786 xclp[2] = xmax; yclp[2] = ymax; 00787 xclp[3] = xmin; yclp[3] = ymax; 00788 xclp[4] = xmin; yclp[4] = ymin; 00789 ( *draw )( xclp, yclp, 5 ); 00790 00791 if ( xclp != _xclp ) 00792 { 00793 free( xclp ); 00794 free( yclp ); 00795 } 00796 00797 return; 00798 } 00799 } 00800 00801 /* Now handle cases where fill polygon intersects two sides of the box */ 00802 00803 if ( iclp >= 2 ) 00804 { 00805 int debug = 0; 00806 int dir = circulation( x, y, npts ); 00807 if ( debug ) 00808 { 00809 if ( ( xclp[0] == xmin && xclp[iclp - 1] == xmax ) || 00810 ( xclp[0] == xmax && xclp[iclp - 1] == xmin ) || 00811 ( yclp[0] == ymin && yclp[iclp - 1] == ymax ) || 00812 ( yclp[0] == ymax && yclp[iclp - 1] == ymin ) || 00813 ( xclp[0] == xmin && yclp[iclp - 1] == ymin ) || 00814 ( yclp[0] == ymin && xclp[iclp - 1] == xmin ) || 00815 ( xclp[0] == xmax && yclp[iclp - 1] == ymin ) || 00816 ( yclp[0] == ymin && xclp[iclp - 1] == xmax ) || 00817 ( xclp[0] == xmax && yclp[iclp - 1] == ymax ) || 00818 ( yclp[0] == ymax && xclp[iclp - 1] == xmax ) || 00819 ( xclp[0] == xmin && yclp[iclp - 1] == ymax ) || 00820 ( yclp[0] == ymax && xclp[iclp - 1] == xmin ) ) 00821 { 00822 printf( "dir=%d, clipped points:\n", dir ); 00823 for ( i = 0; i < iclp; i++ ) 00824 printf( " x[%d]=%d y[%d]=%d", i, xclp[i], i, yclp[i] ); 00825 printf( "\n" ); 00826 printf( "pre-clipped points:\n" ); 00827 for ( i = 0; i < npts; i++ ) 00828 printf( " x[%d]=%d y[%d]=%d", i, x[i], i, y[i] ); 00829 printf( "\n" ); 00830 } 00831 } 00832 00833 /* The cases where the fill region is divided 2/2 */ 00834 /* Divided horizontally */ 00835 if ( xclp[0] == xmin && xclp[iclp - 1] == xmax ) 00836 { 00837 if ( dir > 0 ) 00838 { 00839 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00840 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00841 } 00842 else 00843 { 00844 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00845 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00846 } 00847 } 00848 else if ( xclp[0] == xmax && xclp[iclp - 1] == xmin ) 00849 { 00850 if ( dir > 0 ) 00851 { 00852 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00853 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00854 } 00855 else 00856 { 00857 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00858 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00859 } 00860 } 00861 00862 /* Divided vertically */ 00863 else if ( yclp[0] == ymin && yclp[iclp - 1] == ymax ) 00864 { 00865 if ( dir > 0 ) 00866 { 00867 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00868 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00869 } 00870 else 00871 { 00872 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00873 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00874 } 00875 } 00876 else if ( yclp[0] == ymax && yclp[iclp - 1] == ymin ) 00877 { 00878 if ( dir > 0 ) 00879 { 00880 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00881 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00882 } 00883 else 00884 { 00885 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00886 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00887 } 00888 } 00889 00890 /* The cases where the fill region is divided 3/1 -- 00891 * LL LR UR UL 00892 +-----+ +-----+ +-----+ +-----+ 00893 | | | | | \| |/ | 00894 | | | | | | | | 00895 |\ | | /| | | | | 00896 +-----+ +-----+ +-----+ +-----+ 00897 * 00898 * Note when we go the long way around, if the direction is reversed the 00899 * three vertices must be visited in the opposite order. 00900 */ 00901 /* LL, short way around */ 00902 else if ( ( xclp[0] == xmin && yclp[iclp - 1] == ymin && dir < 0 ) || 00903 ( yclp[0] == ymin && xclp[iclp - 1] == xmin && dir > 0 ) ) 00904 { 00905 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00906 } 00907 /* LL, long way around, counterclockwise */ 00908 else if ( ( xclp[0] == xmin && yclp[iclp - 1] == ymin && dir > 0 ) ) 00909 { 00910 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00911 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00912 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00913 } 00914 /* LL, long way around, clockwise */ 00915 else if ( ( yclp[0] == ymin && xclp[iclp - 1] == xmin && dir < 0 ) ) 00916 { 00917 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00918 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00919 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00920 } 00921 /* LR, short way around */ 00922 else if ( ( xclp[0] == xmax && yclp[iclp - 1] == ymin && dir > 0 ) || 00923 ( yclp[0] == ymin && xclp[iclp - 1] == xmax && dir < 0 ) ) 00924 { 00925 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00926 } 00927 /* LR, long way around, counterclockwise */ 00928 else if ( yclp[0] == ymin && xclp[iclp - 1] == xmax && dir > 0 ) 00929 { 00930 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00931 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00932 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00933 } 00934 /* LR, long way around, clockwise */ 00935 else if ( xclp[0] == xmax && yclp[iclp - 1] == ymin && dir < 0 ) 00936 { 00937 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00938 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00939 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00940 } 00941 /* UR, short way around */ 00942 else if ( ( xclp[0] == xmax && yclp[iclp - 1] == ymax && dir < 0 ) || 00943 ( yclp[0] == ymax && xclp[iclp - 1] == xmax && dir > 0 ) ) 00944 { 00945 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00946 } 00947 /* UR, long way around, counterclockwise */ 00948 else if ( xclp[0] == xmax && yclp[iclp - 1] == ymax && dir > 0 ) 00949 { 00950 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00951 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00952 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00953 } 00954 /* UR, long way around, clockwise */ 00955 else if ( yclp[0] == ymax && xclp[iclp - 1] == xmax && dir < 0 ) 00956 { 00957 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00958 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00959 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00960 } 00961 /* UL, short way around */ 00962 else if ( ( xclp[0] == xmin && yclp[iclp - 1] == ymax && dir > 0 ) || 00963 ( yclp[0] == ymax && xclp[iclp - 1] == xmin && dir < 0 ) ) 00964 { 00965 xclp[iclp] = xmin; yclp[iclp] = ymax; iclp++; 00966 } 00967 /* UL, long way around, counterclockwise */ 00968 else if ( yclp[0] == ymax && xclp[iclp - 1] == xmin && dir > 0 ) 00969 { 00970 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00971 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00972 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00973 } 00974 /* UL, long way around, clockwise */ 00975 else if ( xclp[0] == xmin && yclp[iclp - 1] == ymax && dir < 0 ) 00976 { 00977 xclp[iclp] = xmax; yclp[iclp] = ymax; iclp++; 00978 xclp[iclp] = xmax; yclp[iclp] = ymin; iclp++; 00979 xclp[iclp] = xmin; yclp[iclp] = ymin; iclp++; 00980 } 00981 } 00982 00983 /* Check for the case that only one side has been crossed 00984 * (AM) Just checking a single point turns out not to be 00985 * enough, apparently the crossed_*1 and crossed_*2 variables 00986 * are not quite what I expected. 00987 */ 00988 if ( inside_lb + inside_rb + inside_lu + inside_ru == 4 ) 00989 { 00990 int dir = circulation( x, y, npts ); 00991 PLINT xlim[4], ylim[4]; 00992 int insert; 00993 int incr; 00994 00995 xlim[0] = xmin; ylim[0] = ymin; 00996 xlim[1] = xmax; ylim[1] = ymin; 00997 xlim[2] = xmax; ylim[2] = ymax; 00998 xlim[3] = xmin; ylim[3] = ymax; 00999 01000 if ( crossed_left + crossed_right + crossed_down + crossed_up == 1 ) 01001 { 01002 if ( dir > 0 ) 01003 { 01004 incr = 1; 01005 insert = 0 * crossed_left + 1 * crossed_down + 2 * crossed_right + 01006 3 * crossed_up; 01007 } 01008 else 01009 { 01010 incr = -1; 01011 insert = 3 * crossed_left + 2 * crossed_up + 1 * crossed_right + 01012 0 * crossed_down; 01013 } 01014 } 01015 01016 if ( crossed_left + crossed_right == 2 && crossed_down + crossed_up == 0 ) 01017 { 01018 if ( xclp[iclp - 1] == xmin ) 01019 { 01020 if ( dir == 1 ) 01021 { 01022 incr = 1; 01023 insert = 0; 01024 } 01025 else 01026 { 01027 incr = -1; 01028 insert = 3; 01029 } 01030 } 01031 else 01032 { 01033 if ( dir == 1 ) 01034 { 01035 incr = 1; 01036 insert = 1; 01037 } 01038 else 01039 { 01040 incr = -1; 01041 insert = 2; 01042 } 01043 } 01044 } 01045 01046 if ( crossed_left + crossed_right == 0 && crossed_down + crossed_up == 2 ) 01047 { 01048 if ( yclp[iclp - 1] == ymin ) 01049 { 01050 if ( dir == 1 ) 01051 { 01052 incr = 1; 01053 insert = 1; 01054 } 01055 else 01056 { 01057 incr = -1; 01058 insert = 0; 01059 } 01060 } 01061 else 01062 { 01063 if ( dir == 1 ) 01064 { 01065 incr = 1; 01066 insert = 3; 01067 } 01068 else 01069 { 01070 incr = -1; 01071 insert = 2; 01072 } 01073 } 01074 } 01075 01076 for ( i = 0; i < 4; i++ ) 01077 { 01078 xclp[iclp] = xlim[insert]; 01079 yclp[iclp] = ylim[insert]; 01080 iclp++; 01081 insert += incr; 01082 if ( insert > 3 ) 01083 insert = 0; 01084 if ( insert < 0 ) 01085 insert = 3; 01086 } 01087 } 01088 01089 /* Draw the sucker */ 01090 if ( iclp >= 3 ) 01091 ( *draw )( xclp, yclp, iclp ); 01092 01093 if ( xclp != _xclp ) 01094 { 01095 free( xclp ); 01096 free( yclp ); 01097 } 01098 } 01099 #endif /* USE_FILL_INTERSECTION_POLYGON */ 01100 01101 /*----------------------------------------------------------------------*\ 01102 * int circulation() 01103 * 01104 * Returns the circulation direction for a given polyline: positive is 01105 * counterclockwise, negative is clockwise (right hand rule). 01106 * 01107 * Used to get the circulation of the fill polygon around the bounding box, 01108 * when the fill polygon is larger than the bounding box. Counts left 01109 * (positive) vs right (negative) hand turns using a cross product, instead of 01110 * performing all the expensive trig calculations needed to get this 100% 01111 * correct. For the fill cases encountered in plplot, this treatment should 01112 * give the correct answer most of the time, by far. When used with plshades, 01113 * the typical return value is 3 or -3, since 3 turns are necessary in order 01114 * to complete the fill region. Only for really oddly shaped fill regions 01115 * will it give the wrong answer. 01116 * 01117 * AM: 01118 * Changed the computation: use the outer product to compute the surface 01119 * area, the sign determines if the polygon is followed clockwise or 01120 * counterclockwise. This is more reliable. Floating-point numbers 01121 * are used to avoid overflow. 01122 \*----------------------------------------------------------------------*/ 01123 01124 int 01125 circulation( PLINT *x, PLINT *y, PLINT npts ) 01126 { 01127 PLFLT xproduct; 01128 int direction = 0; 01129 PLFLT x1, y1, x2, y2, x3, y3; 01130 int i; 01131 01132 xproduct = 0.0; 01133 x1 = x[0]; 01134 y1 = y[0]; 01135 for ( i = 1; i < npts - 2; i++ ) 01136 { 01137 x2 = x[i + 1]; 01138 y2 = y[i + 1]; 01139 x3 = x[i + 2]; 01140 y3 = y[i + 2]; 01141 xproduct = xproduct + ( x2 - x1 ) * ( y3 - y2 ) - ( y2 - y1 ) * ( x3 - x2 ); 01142 } 01143 01144 if ( xproduct > 0.0 ) 01145 direction = 1; 01146 if ( xproduct < 0.0 ) 01147 direction = -1; 01148 return direction; 01149 } 01150 01151 01152 /* PLFLT wrapper for !notpointinpolygon. */ 01153 int 01154 plP_pointinpolygon( PLINT n, const PLFLT *x, const PLFLT *y, PLFLT xp, PLFLT yp ) 01155 { 01156 int i, return_value; 01157 PLINT *xint, *yint; 01158 PLFLT xmaximum = fabs( xp ), ymaximum = fabs( yp ), xscale, yscale; 01159 if ( ( xint = (PLINT *) malloc( n * sizeof ( PLINT ) ) ) == NULL ) 01160 { 01161 plexit( "PlP_pointinpolygon: Insufficient memory" ); 01162 } 01163 if ( ( yint = (PLINT *) malloc( n * sizeof ( PLINT ) ) ) == NULL ) 01164 { 01165 plexit( "PlP_pointinpolygon: Insufficient memory" ); 01166 } 01167 for ( i = 0; i < n; i++ ) 01168 { 01169 xmaximum = MAX( xmaximum, fabs( x[i] ) ); 01170 ymaximum = MAX( ymaximum, fabs( y[i] ) ); 01171 } 01172 xscale = 1.e8 / xmaximum; 01173 yscale = 1.e8 / ymaximum; 01174 for ( i = 0; i < n; i++ ) 01175 { 01176 xint[i] = (PLINT) ( xscale * x[i] ); 01177 yint[i] = (PLINT) ( yscale * y[i] ); 01178 } 01179 return_value = !notpointinpolygon( n, xint, yint, 01180 (PLINT) ( xscale * xp ), (PLINT) ( yscale * yp ) ); 01181 free( xint ); 01182 free( yint ); 01183 return return_value; 01184 } 01185 /*----------------------------------------------------------------------*\ 01186 * int notpointinpolygon() 01187 * 01188 * Returns 0, 1, or 2 depending on whether the test point is definitely 01189 * inside, near the border, or definitely outside the polygon. 01190 * Notes: 01191 * This "Ray casting algorithm" has been described in 01192 * http://en.wikipedia.org/wiki/Point_in_polygon. 01193 * Logic still needs to be inserted to take care of the "ray passes 01194 * through vertex" problem in a numerically robust way. 01195 \*----------------------------------------------------------------------*/ 01196 01197 /* Temporary until get rid of old code altogether. */ 01198 #define NEW_NOTPOINTINPOLYGON_CODE 01199 static int 01200 notpointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) 01201 { 01202 #ifdef NEW_NOTPOINTINPOLYGON_CODE 01203 int i, im1, ifnotcrossed; 01204 int count_crossings = 0; 01205 PLINT xmin, xout, yout, xintersect, yintersect; 01206 01207 01208 /* Determine a point outside the polygon */ 01209 01210 xmin = x[0]; 01211 xout = x[0]; 01212 yout = y[0]; 01213 for ( i = 1; i < n; i++ ) 01214 { 01215 xout = MAX( xout, x[i] ); 01216 xmin = MIN( xmin, x[i] ); 01217 } 01218 /* + 10 to make sure completely outside. */ 01219 xout = xout + ( xout - xmin ) + 10; 01220 01221 /* Determine whether the line between (xout, yout) and (xp, yp) intersects 01222 * one of the polygon segments. */ 01223 01224 im1 = n - 1; 01225 for ( i = 0; i < n; i++ ) 01226 { 01227 if ( !( x[im1] == x[i] && y[im1] == y[i] ) ) 01228 { 01229 ifnotcrossed = notcrossed( &xintersect, &yintersect, 01230 x[im1], y[im1], x[i], y[i], 01231 xp, yp, xout, yout ); 01232 01233 if ( !ifnotcrossed ) 01234 count_crossings++; 01235 else if ( ifnotcrossed & ( PL_NEAR_A1 | PL_NEAR_A2 | PL_NEAR_B1 | PL_NEAR_B2 ) ) 01236 return 1; 01237 } 01238 im1 = i; 01239 } 01240 01241 /* return 0 if the test point is definitely inside 01242 * (count_crossings odd), return 1 if the test point is near (see 01243 * above logic), and return 2 if the test point is definitely 01244 * outside the border (count_crossings even). */ 01245 if ( ( count_crossings % 2 ) == 1 ) 01246 return 0; 01247 else 01248 return 2; 01249 } 01250 #else /* NEW_NOTPOINTINPOLYGON_CODE */ 01251 int i; 01252 int count_crossings; 01253 PLFLT x1, y1, x2, y2, xpp, ypp, xout, yout, xmax; 01254 PLFLT xvp, yvp, xvv, yvv, xv1, yv1, xv2, yv2; 01255 PLFLT inprod1, inprod2; 01256 01257 xpp = (PLFLT) xp; 01258 ypp = (PLFLT) yp; 01259 01260 count_crossings = 0; 01261 01262 01263 /* Determine a point outside the polygon */ 01264 01265 xmax = x[0]; 01266 xout = x[0]; 01267 yout = y[0]; 01268 for ( i = 0; i < n; i++ ) 01269 { 01270 if ( xout > x[i] ) 01271 { 01272 xout = x[i]; 01273 } 01274 if ( xmax < x[i] ) 01275 { 01276 xmax = x[i]; 01277 } 01278 } 01279 xout = xout - ( xmax - xout ); 01280 01281 /* Determine for each side whether the line segment between 01282 * our two points crosses the vertex */ 01283 01284 xpp = (PLFLT) xp; 01285 ypp = (PLFLT) yp; 01286 01287 xvp = xpp - xout; 01288 yvp = ypp - yout; 01289 01290 for ( i = 0; i < n; i++ ) 01291 { 01292 x1 = (PLFLT) x[i]; 01293 y1 = (PLFLT) y[i]; 01294 if ( i < n - 1 ) 01295 { 01296 x2 = (PLFLT) x[i + 1]; 01297 y2 = (PLFLT) y[i + 1]; 01298 } 01299 else 01300 { 01301 x2 = (PLFLT) x[0]; 01302 y2 = (PLFLT) y[0]; 01303 } 01304 01305 /* Skip zero-length segments */ 01306 if ( x1 == x2 && y1 == y2 ) 01307 { 01308 continue; 01309 } 01310 01311 /* Line through the two fixed points: 01312 * Are x1 and x2 on either side? */ 01313 xv1 = x1 - xout; 01314 yv1 = y1 - yout; 01315 xv2 = x2 - xout; 01316 yv2 = y2 - yout; 01317 inprod1 = xv1 * yvp - yv1 * xvp; /* Well, with the normal vector */ 01318 inprod2 = xv2 * yvp - yv2 * xvp; 01319 if ( inprod1 * inprod2 >= 0.0 ) 01320 { 01321 /* No crossing possible! */ 01322 continue; 01323 } 01324 01325 /* Line through the two vertices: 01326 * Are xout and xpp on either side? */ 01327 xvv = x2 - x1; 01328 yvv = y2 - y1; 01329 xv1 = xpp - x1; 01330 yv1 = ypp - y1; 01331 xv2 = xout - x1; 01332 yv2 = yout - y1; 01333 inprod1 = xv1 * yvv - yv1 * xvv; 01334 inprod2 = xv2 * yvv - yv2 * xvv; 01335 if ( inprod1 * inprod2 >= 0.0 ) 01336 { 01337 /* No crossing possible! */ 01338 continue; 01339 } 01340 01341 /* We do have a crossing */ 01342 count_crossings++; 01343 } 01344 01345 /* Return the result: an even number of crossings means the 01346 * point is outside the polygon */ 01347 01348 return !( count_crossings % 2 ); 01349 } 01350 #endif /* NEW_NOTPOINTINPOLYGON_CODE */ 01351 01352 #define MAX_RECURSION_DEPTH 10 01353 01354 /* Fill intersection of two simple polygons that do no self-intersect, 01355 * and which have no duplicate vertices or two consecutive edges that 01356 * are parallel. A further requirement is that both polygons have a 01357 * positive orientation (see 01358 * http://en.wikipedia.org/wiki/Curve_orientation). That is, as you 01359 * traverse the boundary in index order, the inside area of the 01360 * polygon is always on the left. Finally, the first vertex of 01361 * polygon 1 (starting with n1 -1) that is not near the border of 01362 * polygon 2 must be outside polygon 2. N.B. it is the calling 01363 * routine's responsibility to insure all those requirements are 01364 * satisfied. 01365 * 01366 * Two polygons that do not self intersect must have an even number of 01367 * edge crossings between them. (ignoring vertex intersections which 01368 * touch, but do not cross). fill_intersection_polygon eliminates 01369 * those intersection crossings by recursion (calling the same routine 01370 * twice again with the second polygon split at a boundary defined by 01371 * the first intersection point, all polygon 1 vertices between the 01372 * intersections, and the second intersection point). Once the 01373 * recursion has eliminated all crossing edges, fill or not using the 01374 * appropriate polygon depending on whether the first and second 01375 * polygons are identical or whether one of them is entirely inside 01376 * the other of them. If ifextrapolygon is true, the fill step will 01377 * consist of another recursive call to the routine with 01378 * ifextrapolygon false, and the second polygon set to an additional 01379 * polygon defined by the stream (not yet implemented). */ 01380 01381 /* arguments to intersection_polygon: 01382 * recursion_depth is just what it says. 01383 * ifextrapolygon used to decide whether to use extra polygon from the stream. 01384 *fill is the fill routine. 01385 *x1, *y1, n1 define the polygon 1 vertices. 01386 * i1start is the first polygon 1 index to look at (because all the previous 01387 * ones have been completely processed). 01388 *x2, *y2, *if2, n2 define the polygon 2 vertices plus a status indicator 01389 * for each vertex which is 1 for a previous crossing and 2 for a polygon 01390 * 1 vertex. 01391 * fill_status is 1 when polygons 1 and 2 _must_ include some joint 01392 * filled area and is -1 when polygons 1 and 2 _must_ include some 01393 * unfilled area. fill_status of +/- 1 is determined from the 01394 * orientations of polygon 1 and 2 from the next higher recursion 01395 * level and how those two are combined to form the polygon 2 01396 * split at this recursion level. fill_status = 0 occurs (at 01397 * recursion level 0) for polygons 1 and 2 that are independent of 01398 * each other. */ 01399 01400 void 01401 fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, 01402 PLINT fill_status, 01403 void ( *fill )( short *, short *, PLINT ), 01404 const PLINT *x1, const PLINT *y1, 01405 PLINT i1start, PLINT n1, 01406 const PLINT *x2, const PLINT *y2, 01407 const PLINT *if2, PLINT n2 ) 01408 { 01409 PLINT i1, i1m1, i1start_new, 01410 i2, i2m1, 01411 kk, kkstart1, kkstart21, kkstart22, 01412 k, kstart, range1, 01413 range21, range22, ncrossed, ncrossed_change, 01414 nsplit1, nsplit2, nsplit2m1; 01415 PLINT xintersect[2], yintersect[2], i1intersect[2], 01416 i2intersect[2], ifnotcrossed; 01417 PLINT *xsplit1, *ysplit1, *ifsplit1, 01418 *xsplit2, *ysplit2, *ifsplit2; 01419 PLINT ifill, nfill = 0, 01420 ifnotpolygon1inpolygon2, ifnotpolygon2inpolygon1; 01421 const PLINT *xfiller, *yfiller; 01422 short *xfill, *yfill; 01423 01424 if ( recursion_depth > MAX_RECURSION_DEPTH ) 01425 { 01426 plwarn( "fill_intersection_polygon: Recursion_depth too large. " 01427 "Probably an internal error figuring out intersections. " ); 01428 return; 01429 } 01430 01431 if ( n1 < 3 ) 01432 { 01433 plwarn( "fill_intersection_polygon: Internal error; n1 < 3." ); 01434 return; 01435 } 01436 01437 if ( n2 < 3 ) 01438 { 01439 plwarn( "fill_intersection_polygon: Internal error; n2 < 3." ); 01440 return; 01441 } 01442 01443 if ( i1start < 0 || i1start >= n1 ) 01444 { 01445 plwarn( "fill_intersection_polygon: invalid i1start." ); 01446 return; 01447 } 01448 01449 /* Check that there are no duplicate vertices. */ 01450 i1m1 = i1start - 1; 01451 if ( i1m1 < 0 ) 01452 i1m1 = n1 - 1; 01453 01454 for ( i1 = i1start; i1 < n1; i1++ ) 01455 { 01456 if ( x1[i1] == x1[i1m1] && y1[i1] == y1[i1m1] ) 01457 break; 01458 i1m1 = i1; 01459 } 01460 01461 if ( i1 < n1 ) 01462 { 01463 plwarn( "fill_intersection_polygon: Internal error; i1 < n1." ); 01464 return; 01465 } 01466 01467 i2m1 = n2 - 1; 01468 for ( i2 = 0; i2 < n2; i2++ ) 01469 { 01470 if ( x2[i2] == x2[i2m1] && y2[i2] == y2[i2m1] ) 01471 break; 01472 i2m1 = i2; 01473 } 01474 01475 if ( i2 < n2 ) 01476 { 01477 plwarn( "fill_intersection_polygon: Internal error; i2 < n2." ); 01478 return; 01479 } 01480 01481 /* 01482 * 01483 * Follow polygon 1 (checking intersections with polygon 2 for each 01484 * segment of polygon 1) until you have accumulated two 01485 * intersections with polygon 2. Here is an ascii-art illustration 01486 * of the situation. 01487 * 01488 * 01489 * 2???2 01490 * 01491 * 2 2 01492 * 01493 * --- 1 1 01494 * 1 2 1 1 ... 01495 * X 01496 * 1 01497 * X 01498 * 2 01499 * 1 1 01500 * 1 01501 * 2 01502 * 2 01503 * 2???2 01504 * 01505 * 01506 * "1" marks polygon 1 vertices, "2" marks polygon 2 vertices, "X" 01507 * marks the intersections, "---" stands for part of polygon 1 01508 * that has been previously searched for all possible intersections 01509 * from index 0, and "..." means polygon 1 continues 01510 * with more potential intersections above and/or below this diagram 01511 * before it finally hooks back to connect with the index 0 vertex. 01512 * "2???2" stands for parts of polygon 2 that must connect with each other 01513 * (since the polygon 1 path between the two intersections is 01514 * known to be free of intersections.) 01515 * 01516 * Polygon 2 is split at the boundary defined by the two 01517 * intersections and all (in this case three) polygon 1 vertices 01518 * between the two intersections for the next recursion level. We 01519 * absolutely know for that boundary that no more intersections can 01520 * occur (both polygon 1 and polygon 2 are guaranteed not to 01521 * self-intersect) so we mark the status of those vertices with that 01522 * information so those polygon 2 split vertices will not be used to 01523 * search for further intersections at deeper recursion levels. 01524 * Note, we know nothing about whether the remaining "2???2" parts of the 01525 * split polygon 2 intersect with polygon 1 or not so those will 01526 * continued to be searched at deeper recursion levels. At the same 01527 * time, we absolutely know that the part of polygon 1 to the left of 01528 * rightmost x down to and including index 0 cannot yield more 01529 * intersections with any split of polygon 2 so we adjust the lower 01530 * limit of polygon 1 to be used for intersection searches at deeper 01531 * recursion levels. The result should be that at sufficiently deep 01532 * recursion depth we always end up with the case that there are no 01533 * intersections to be found between polygon 1 and some polygon 2 01534 * split, and in that case we move on to the end phase below. 01535 */ 01536 ncrossed = 0; 01537 i1m1 = i1start - 1; 01538 if ( i1m1 < 0 ) 01539 i1m1 += n1; 01540 for ( i1 = i1start; i1 < n1; i1++ ) 01541 { 01542 ncrossed_change = number_crossings( 01543 &xintersect[ncrossed], &yintersect[ncrossed], 01544 &i2intersect[ncrossed], 2 - ncrossed, 01545 i1, n1, x1, y1, n2, x2, y2 ); 01546 if ( ncrossed_change > 0 ) 01547 { 01548 i1intersect[ncrossed] = i1; 01549 if ( ncrossed_change == 2 ) 01550 ; 01551 i1intersect[1] = i1; 01552 01553 ncrossed += ncrossed_change; 01554 if ( ncrossed == 2 ) 01555 { 01556 /* Have discovered the first two crossings for 01557 * polygon 1 at i1 = i1start or above. */ 01558 01559 /* New i1start is the same as the current i1 (just 01560 * in case there are more crossings to find between 01561 * i1m1 and i1.) */ 01562 i1start_new = i1; 01563 01564 /* Split polygon 2 at the boundary consisting of 01565 * first intersection, intervening (if any) range1 01566 * polygon 1 points and second intersection. */ 01567 /* range1 must always be non-negative because i1 01568 * range only traversed once. */ 01569 range1 = i1intersect[1] - i1intersect[0]; 01570 /* Polygon 2 intersects could be anywhere (since 01571 * i2 range repeated until get an intersect). 01572 * Divide polygon 2 into two polygons with a 01573 * common boundary consisting of the first intersect, 01574 * range1 points from polygon 1 starting at index 01575 * kkstart1 of polygon 1, and the second intersect. */ 01576 kkstart1 = i1intersect[0]; 01577 01578 /* Calculate polygon 2 index range in split 1 (the 01579 * split that proceeds beyond the second intersect with 01580 * ascending i2 values). */ 01581 range21 = i2intersect[0] - i2intersect[1]; 01582 if ( range21 < 0 ) 01583 range21 += n2; 01584 /* i2 intersect values range between 0 and n2 - 1 so 01585 * the smallest untransformed range21 value is -n2 + 1, 01586 * and the largest untransformed range21 value is n2 - 1. 01587 * This means the smallest transformed range21 value is 0 01588 * (which occurs only ifi2intersect[0] = i2intersect[1], 01589 * see more commentary for that special case below) while 01590 * the largest transformed range21 value is n2 - 1. */ 01591 01592 if ( range21 == 0 ) 01593 { 01594 int ifxsort, ifascend; 01595 /* For this case, the two crossings occur within the same 01596 * polygon 2 boundary segment and if those two crossings 01597 * are in ascending/descending order in i2, then split 1 01598 * (the split with the positive fill_status) must include 01599 * all/none of the points in polygon 2. */ 01600 i2 = i2intersect[1]; 01601 i2m1 = i2 - 1; 01602 if ( i2m1 < 0 ) 01603 i2m1 += n2; 01604 01605 ifxsort = abs( x2[i2] - x2[i2m1] ) > abs( y2[i2] - y2[i2m1] ); 01606 ifascend = ( ifxsort && x2[i2] > x2[i2m1] ) || 01607 ( !ifxsort && y2[i2] > y2[i2m1] ); 01608 if ( ( ifxsort && ifascend && xintersect[0] < xintersect[1] ) || 01609 ( !ifxsort && ifascend && yintersect[0] < yintersect[1] ) || 01610 ( ifxsort && !ifascend && xintersect[0] >= xintersect[1] ) || 01611 ( !ifxsort && !ifascend && yintersect[0] >= yintersect[1] ) ) 01612 { 01613 range21 = n2; 01614 } 01615 } 01616 01617 kkstart21 = i2intersect[1]; 01618 nsplit1 = 2 + range1 + range21; 01619 01620 /* Split 2 of polygon 2 consists of the 01621 * boundary + range22 (= n2 - range21) points 01622 * between kkstart22 (= i2intersect[1]-1) and i2intersect[0] in 01623 * descending order of polygon 2 indices. */ 01624 range22 = n2 - range21; 01625 /* Starting i2 index of split 2. */ 01626 kkstart22 = i2intersect[1] - 1; 01627 if ( kkstart22 < 0 ) 01628 kkstart22 += n2; 01629 nsplit2 = 2 + range1 + range22; 01630 01631 if ( ( xsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ) ) ) == NULL ) 01632 { 01633 plexit( "fill_intersection_polygon: Insufficient memory" ); 01634 } 01635 if ( ( ysplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ) ) ) == NULL ) 01636 { 01637 plexit( "fill_intersection_polygon: Insufficient memory" ); 01638 } 01639 if ( ( ifsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ) ) ) == NULL ) 01640 { 01641 plexit( "fill_intersection_polygon: Insufficient memory" ); 01642 } 01643 01644 if ( ( xsplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ) ) ) == NULL ) 01645 { 01646 plexit( "fill_intersection_polygon: Insufficient memory" ); 01647 } 01648 if ( ( ysplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ) ) ) == NULL ) 01649 { 01650 plexit( "fill_intersection_polygon: Insufficient memory" ); 01651 } 01652 if ( ( ifsplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ) ) ) == NULL ) 01653 { 01654 plexit( "fill_intersection_polygon: Insufficient memory" ); 01655 } 01656 /* Common boundary between split1 and split2. */ 01657 /* N.B. Although basic index arithmetic for 01658 * split 2 is done in negative orientation 01659 * order because the index is decrementing 01660 * relative to the index of split 2, actually 01661 * store results in reverse order to preserve 01662 * the positive orientation that by assumption 01663 * both polygon 1 and 2 have. */ 01664 k = 0; 01665 xsplit1[k] = xintersect[0]; 01666 ysplit1[k] = yintersect[0]; 01667 ifsplit1[k] = 1; 01668 nsplit2m1 = nsplit2 - 1; 01669 xsplit2[nsplit2m1 - k] = xintersect[0]; 01670 ysplit2[nsplit2m1 - k] = yintersect[0]; 01671 ifsplit2[nsplit2m1 - k] = 1; 01672 kstart = k + 1; 01673 kk = kkstart1; 01674 /* No wrap checks on kk index below because 01675 * it must always be in valid range (since 01676 * polygon 1 traversed only once). */ 01677 for ( k = kstart; k < range1 + 1; k++ ) 01678 { 01679 xsplit1[k] = x1[kk]; 01680 ysplit1[k] = y1[kk]; 01681 ifsplit1[k] = 2; 01682 xsplit2[nsplit2m1 - k] = x1[kk]; 01683 ysplit2[nsplit2m1 - k] = y1[kk++]; 01684 ifsplit2[nsplit2m1 - k] = 2; 01685 } 01686 xsplit1[k] = xintersect[1]; 01687 ysplit1[k] = yintersect[1]; 01688 ifsplit1[k] = 1; 01689 xsplit2[nsplit2m1 - k] = xintersect[1]; 01690 ysplit2[nsplit2m1 - k] = yintersect[1]; 01691 ifsplit2[nsplit2m1 - k] = 1; 01692 01693 /* Finish off collecting split1 using ascending kk 01694 * values. */ 01695 kstart = k + 1; 01696 kk = kkstart21; 01697 for ( k = kstart; k < nsplit1; k++ ) 01698 { 01699 xsplit1[k] = x2[kk]; 01700 ysplit1[k] = y2[kk]; 01701 ifsplit1[k] = if2[kk++]; 01702 if ( kk >= n2 ) 01703 kk -= n2; 01704 } 01705 01706 /* N.B. the positive orientation of split1 is 01707 * preserved since the index order is the same 01708 * as that of polygon 2, and by assumption 01709 * that polygon and polygon 1 have identical 01710 * positive orientations. */ 01711 fill_intersection_polygon( 01712 recursion_depth + 1, ifextrapolygon, 1, fill, 01713 x1, y1, i1start_new, n1, 01714 xsplit1, ysplit1, ifsplit1, nsplit1 ); 01715 free( xsplit1 ); 01716 free( ysplit1 ); 01717 free( ifsplit1 ); 01718 01719 /* Finish off collecting split2 using descending kk 01720 * values. */ 01721 kk = kkstart22; 01722 for ( k = kstart; k < nsplit2; k++ ) 01723 { 01724 xsplit2[nsplit2m1 - k] = x2[kk]; 01725 ysplit2[nsplit2m1 - k] = y2[kk]; 01726 ifsplit2[nsplit2m1 - k] = if2[kk--]; 01727 if ( kk < 0 ) 01728 kk += n2; 01729 } 01730 01731 /* N.B. the positive orientation of split2 is 01732 * preserved since the index order is the same 01733 * as that of polygon 2, and by assumption 01734 * that polygon and polygon 1 have identical 01735 * positive orientations. */ 01736 fill_intersection_polygon( 01737 recursion_depth + 1, ifextrapolygon, -1, fill, 01738 x1, y1, i1start_new, n1, 01739 xsplit2, ysplit2, ifsplit2, nsplit2 ); 01740 free( xsplit2 ); 01741 free( ysplit2 ); 01742 free( ifsplit2 ); 01743 return; 01744 } 01745 } 01746 i1m1 = i1; 01747 } 01748 01749 if ( ncrossed != 0 ) 01750 { 01751 plwarn( "fill_intersection_polygon: Internal error; ncrossed != 0." ); 01752 return; 01753 } 01754 01755 /* This end phase is reached only if no crossings are found. */ 01756 01757 /* If a fill_status of +/- 1 is known, use that to fill or not since 01758 +1 corresponds to all of polygon 2 inside polygon 1 and -1 01759 * corresponds to none of polygon 2 inside polygon 1. */ 01760 if ( fill_status == -1 ) 01761 return; 01762 else if ( fill_status == 1 ) 01763 { 01764 nfill = n2; 01765 xfiller = x2; 01766 yfiller = y2; 01767 } 01768 else if ( fill_status == 0 ) 01769 //else if ( 1 ) 01770 { 01771 if ( recursion_depth != 0 ) 01772 { 01773 plwarn( "fill_intersection_polygon: Internal error; fill_status == 0 for recursion_depth > 0" ); 01774 return; 01775 } 01776 /* For this case (recursion level 0) the two polygons are 01777 * completely independent with no crossings between them or 01778 * edges constructed from one another. 01779 * 01780 * The intersection of polygon 2 and 1, must be either of them (in 01781 * which case fill with the inner one), or neither of them (in 01782 * which case don't fill at all). */ 01783 01784 /* Classify polygon 1 by looking for first vertex in polygon 1 01785 * that is definitely inside or outside polygon 2. */ 01786 for ( i1 = 0; i1 < n1; i1++ ) 01787 { 01788 if ( ( ifnotpolygon1inpolygon2 = 01789 notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] ) ) != 1 ) 01790 break; 01791 } 01792 01793 /* Classify polygon 2 by looking for first vertex in polygon 2 01794 * that is definitely inside or outside polygon 1. */ 01795 ifnotpolygon2inpolygon1 = 1; 01796 for ( i2 = 0; i2 < n2; i2++ ) 01797 { 01798 /* Do not bother checking vertices already known to be on the 01799 * boundary with polygon 1. */ 01800 if ( !if2[i2] && ( ifnotpolygon2inpolygon1 = 01801 notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] ) ) != 1 ) 01802 break; 01803 } 01804 01805 if ( ifnotpolygon2inpolygon1 == 0 && ifnotpolygon1inpolygon2 == 0 ) 01806 plwarn( "fill_intersection_polygon: Internal error; no intersections found but each polygon definitely inside the other!" ); 01807 else if ( ifnotpolygon2inpolygon1 == 2 && ifnotpolygon1inpolygon2 == 2 ) 01808 /* The polygons do not intersect each other so do not fill in this 01809 * case. */ 01810 return; 01811 else if ( ifnotpolygon2inpolygon1 == 0 ) 01812 { 01813 /* Polygon 2 definitely inside polygon 1. */ 01814 nfill = n2; 01815 xfiller = x2; 01816 yfiller = y2; 01817 } 01818 else if ( ifnotpolygon1inpolygon2 == 0 ) 01819 { 01820 /* Polygon 1 definitely inside polygon 2. */ 01821 nfill = n1; 01822 xfiller = x1; 01823 yfiller = y1; 01824 } 01825 else if ( ifnotpolygon2inpolygon1 == 1 && ifnotpolygon1inpolygon2 == 1 ) 01826 { 01827 /* Polygon 2 vertices near polygon 1 border and vice versa which 01828 * implies the polygons are identical. */ 01829 nfill = n2; 01830 xfiller = x2; 01831 yfiller = y2; 01832 } 01833 else 01834 { 01835 /* Polygon 1 inscribed in polygon 2 or vice versa. This is normally 01836 * unlikely for two independent polygons so the implementation is 01837 * ToDo. */ 01838 plwarn( "fill_intersection_polygon: inscribed polygons are still ToDo" ); 01839 } 01840 } 01841 01842 if ( nfill > 0 ) 01843 { 01844 if ( ( xfill = (short *) malloc( nfill * sizeof ( short ) ) ) == NULL ) 01845 { 01846 plexit( "fill_intersection_polygon: Insufficient memory" ); 01847 } 01848 if ( ( yfill = (short *) malloc( nfill * sizeof ( short ) ) ) == NULL ) 01849 { 01850 plexit( "fill_intersection_polygon: Insufficient memory" ); 01851 } 01852 for ( ifill = 0; ifill < nfill; ifill++ ) 01853 { 01854 xfill[ifill] = xfiller[ifill]; 01855 yfill[ifill] = yfiller[ifill]; 01856 } 01857 ( *fill )( xfill, yfill, nfill ); 01858 free( xfill ); 01859 free( yfill ); 01860 } 01861 01862 return; 01863 } 01864 01865 /* Returns a 0 status code if the two line segments A, and B defined 01866 * by their end points (xA1, yA1, xA2, yA2, xB1, yB1, xB2, and yB2) 01867 * definitely (i.e., intersection point is not near the ends of either 01868 * of the line segments) cross each other. Otherwise, return a status 01869 * code specifying the type of non-crossing (i.e., completely 01870 * separate, near one of the ends, parallel). Return the actual 01871 * intersection (or average of the x end points and y end points for 01872 * the parallel, not crossed case) via the argument list pointers to 01873 * xintersect and yintersect. */ 01874 01875 int 01876 notcrossed( PLINT * xintersect, PLINT * yintersect, 01877 PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, 01878 PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ) 01879 { 01880 PLFLT factor, factor_NBCC, fxintersect, fyintersect; 01881 /* These variables are PLFLT not for precision, but to 01882 * avoid integer overflows if they were typed as PLINT. */ 01883 PLFLT xA2A1, yA2A1, xB2B1, yB2B1; 01884 PLFLT xB1A1, yB1A1, xB2A1, yB2A1; 01885 PLINT status = 0; 01886 /* 01887 * Two linear equations to be solved for x and y. 01888 * y = ((x - xA1)*yA2 + (xA2 - x)*yA1)/(xA2 - xA1) 01889 * y = ((x - xB1)*yB2 + (xB2 - x)*yB1)/(xB2 - xB1) 01890 * 01891 * Transform those two equations to coordinate system with origin 01892 * at (xA1, yA1). 01893 * y' = x'*yA2A1/xA2A1 01894 * y' = ((x' - xB1A1)*yB2A1 + (xB2A1 - x')*yB1A1)/xB2B1 01895 * ==> 01896 * x' = -( 01897 * (-xB1A1*yB2A1 + xB2A1*yB1A1)/xB2B1)/ 01898 * (yB2B1/xB2B1 - yA2A1/xA2A1) 01899 * = (xB1A1*yB2A1 - xB2A1*yB1A1)*xA2A1/ 01900 * (xA2A1*yB2B1 - yA2A1*xB2B1) 01901 * 01902 */ 01903 01904 xA2A1 = xA2 - xA1; 01905 yA2A1 = yA2 - yA1; 01906 xB2B1 = xB2 - xB1; 01907 yB2B1 = yB2 - yB1; 01908 01909 factor = xA2A1 * yB2B1 - yA2A1 * xB2B1; 01910 factor_NBCC = PL_NBCC * ( fabs( xA2A1 ) + fabs( yB2B1 ) + fabs( yA2A1 ) + fabs( xB2B1 ) ); 01911 if ( fabs( factor ) <= factor_NBCC ) 01912 { 01913 if ( fabs( factor ) > 0. ) 01914 status = status | PL_NEAR_PARALLEL; 01915 else 01916 status = status | PL_PARALLEL; 01917 /* Choice of intersect is arbitrary in this case. Choose A1, A2, 01918 * B1, or B2 (in that order) if any of them lie inside or near 01919 * the other line segment. Otherwise, choose the average point. */ 01920 if ( ( BETW_NBCC( xA1, xB1, xB2 ) && BETW_NBCC( yA1, yB1, yB2 ) ) ) 01921 { 01922 fxintersect = xA1; 01923 fyintersect = yA1; 01924 } 01925 else if ( ( BETW_NBCC( xA2, xB1, xB2 ) && BETW_NBCC( yA2, yB1, yB2 ) ) ) 01926 { 01927 fxintersect = xA2; 01928 fyintersect = yA2; 01929 } 01930 else if ( ( BETW_NBCC( xB1, xA1, xA2 ) && BETW_NBCC( yB1, yA1, yA2 ) ) ) 01931 { 01932 fxintersect = xB1; 01933 fyintersect = yB1; 01934 } 01935 else if ( ( BETW_NBCC( xB2, xA1, xA2 ) && BETW_NBCC( yB2, yA1, yA2 ) ) ) 01936 { 01937 fxintersect = xB2; 01938 fyintersect = yB2; 01939 } 01940 else 01941 { 01942 fxintersect = 0.25 * ( xA1 + xA2 + xB1 + xB2 ); 01943 fyintersect = 0.25 * ( yA1 + yA2 + yB1 + yB2 ); 01944 } 01945 } 01946 else 01947 { 01948 xB1A1 = xB1 - xA1; 01949 yB1A1 = yB1 - yA1; 01950 xB2A1 = xB2 - xA1; 01951 yB2A1 = yB2 - yA1; 01952 01953 factor = ( xB1A1 * yB2A1 - yB1A1 * xB2A1 ) / factor; 01954 fxintersect = factor * xA2A1 + xA1; 01955 fyintersect = factor * yA2A1 + yA1; 01956 } 01957 /* The "redundant" x and y segment range checks (which include near the 01958 * end points) are needed for the vertical and the horizontal cases. */ 01959 if ( ( BETW_NBCC( fxintersect, xA1, xA2 ) && BETW_NBCC( fyintersect, yA1, yA2 ) ) && 01960 ( BETW_NBCC( fxintersect, xB1, xB2 ) && BETW_NBCC( fyintersect, yB1, yB2 ) ) ) 01961 { 01962 /* The intersect is close (within +/- PL_NBCC) to an end point or 01963 * corresponds to a definite crossing of the two line segments. 01964 * Find out which. */ 01965 if ( fabs( fxintersect - xA1 ) <= PL_NBCC && fabs( fyintersect - yA1 ) <= PL_NBCC ) 01966 status = status | PL_NEAR_A1; 01967 else if ( fabs( fxintersect - xA2 ) <= PL_NBCC && fabs( fyintersect - yA2 ) <= PL_NBCC ) 01968 status = status | PL_NEAR_A2; 01969 else if ( fabs( fxintersect - xB1 ) <= PL_NBCC && fabs( fyintersect - yB1 ) <= PL_NBCC ) 01970 status = status | PL_NEAR_B1; 01971 else if ( fabs( fxintersect - xB2 ) <= PL_NBCC && fabs( fyintersect - yB2 ) <= PL_NBCC ) 01972 status = status | PL_NEAR_B2; 01973 /* N.B. if none of the above conditions hold then status remains at 01974 * zero to signal we have a definite crossing. */ 01975 } 01976 else 01977 status = status | PL_NOT_CROSSED; 01978 *xintersect = fxintersect; 01979 *yintersect = fyintersect; 01980 01981 return status; 01982 } 01983 01984 /* Decide if polygon has a positive orientation or not. 01985 * Note the simple algorithm given in 01986 * http://en.wikipedia.org/wiki/Curve_orientation is incorrect for 01987 * non-convex polygons. Instead, for the general nonintersecting case 01988 * use the polygon area method given by 01989 * http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/ where 01990 * you project each edge of the polygon down to the X axis and take the 01991 * area of the enclosed trapezoid. The trapezoid areas outside the 01992 * polygon are completely cancelled if you sum over all edges. Furthermore, 01993 * the sum of the trapezoid areas have terms which are zero when calculated 01994 * with the telescoping rule so the final formula is quite simple. */ 01995 int 01996 positive_orientation( PLINT n, const PLINT *x, const PLINT *y ) 01997 { 01998 PLINT i, im1; 01999 /* Use PLFLT for all calculations to avoid integer overflows. Also, 02000 * the normal PLFLT has 64-bits which means you get exact integer 02001 * arithmetic well beyond the normal integer overflow limits. */ 02002 PLFLT twice_area = 0.; 02003 if ( n < 3 ) 02004 { 02005 plwarn( "positive_orientation: internal logic error, n < 3" ); 02006 return 0; 02007 } 02008 im1 = n - 1; 02009 for ( i = 0; i < n; i++ ) 02010 { 02011 twice_area += (PLFLT) x[im1] * (PLFLT) y[i] - (PLFLT) x[i] * (PLFLT) y[im1]; 02012 im1 = i; 02013 } 02014 if ( twice_area == 0. ) 02015 { 02016 plwarn( "positive_orientation: internal logic error, twice_area == 0." ); 02017 return 0; 02018 } 02019 else 02020 return twice_area > 0.; 02021 } 02022 02023 /* For the line with endpoints which are the (i1-1)th, and i1th 02024 * vertices of polygon 1 with polygon 2 find all definite crossings 02025 * with polygon 1. (The full polygon 1 information is needed only to 02026 * help sort out (NOT IMPLEMENTED YET) ambiguous crossings near a 02027 * vertex of polygon 1.) Sort those definite crossings in ascending 02028 * order along the line between the (i1-1)th and i1th vertices of 02029 * polygon 1, and return the first ncross (= 1 or = 2) crossings in the 02030 * xcross, ycross, and i2cross arrays. Return a zero or positive 02031 * status code of the actual number of crossings retained up to the 02032 * maximum allowed value of ncross. If some error occurred, return a 02033 * negative status code. */ 02034 02035 int 02036 number_crossings( PLINT *xcross, PLINT *ycross, PLINT *i2cross, PLINT ncross, 02037 PLINT i1, PLINT n1, const PLINT *x1, const PLINT *y1, 02038 PLINT n2, const PLINT *x2, const PLINT *y2 ) 02039 { 02040 int i1m1, i2, i2m1, ifnotcrossed; 02041 int ifxsort, ifascend, count_crossings = 0, status = 0; 02042 PLINT xintersect, yintersect; 02043 02044 i1m1 = i1 - 1; 02045 if ( i1m1 < 0 ) 02046 i1m1 += n1; 02047 if ( !( ncross == 1 || ncross == 2 ) || 02048 ( x1[i1m1] == x1[i1] && y1[i1m1] == y1[i1] ) || n1 < 2 || n2 < 2 ) 02049 { 02050 plwarn( "findcrossings: invalid call" ); 02051 return -1; 02052 } 02053 02054 ifxsort = abs( x1[i1] - x1[i1m1] ) > abs( y1[i1] - y1[i1m1] ); 02055 ifascend = ( ifxsort && x1[i1] > x1[i1m1] ) || 02056 ( !ifxsort && y1[i1] > y1[i1m1] ); 02057 02058 /* Determine all crossings between the line between the (i1-1)th 02059 * and i1th vertices of polygon 1 and all edges of polygon 2. 02060 * Retain the lowest and (if ncross ==2) next lowest crossings in 02061 * order of x (or y if ifxsort is false) along the line from i1-1 02062 * to i1. */ 02063 02064 i1m1 = i1 - 1; 02065 if ( i1m1 < 0 ) 02066 i1m1 += n1; 02067 i2m1 = n2 - 1; 02068 for ( i2 = 0; i2 < n2; i2++ ) 02069 { 02070 if ( !( x2[i2m1] == x2[i2] && y2[i2m1] == y2[i2] ) ) 02071 { 02072 ifnotcrossed = notcrossed( &xintersect, &yintersect, 02073 x1[i1m1], y1[i1m1], x1[i1], y1[i1], 02074 x2[i2m1], y2[i2m1], x2[i2], y2[i2] ); 02075 02076 if ( !ifnotcrossed ) 02077 { 02078 count_crossings++; 02079 if ( count_crossings == 1 ) 02080 { 02081 xcross[0] = xintersect; 02082 ycross[0] = yintersect; 02083 i2cross[0] = i2; 02084 status = 1; 02085 } 02086 else 02087 { 02088 if ( ( ifxsort && ifascend && xintersect < xcross[0] ) || 02089 ( !ifxsort && ifascend && yintersect < ycross[0] ) || 02090 ( ifxsort && !ifascend && xintersect >= xcross[0] ) || 02091 ( !ifxsort && !ifascend && yintersect >= ycross[0] ) ) 02092 { 02093 if ( ncross == 2 ) 02094 { 02095 xcross[1] = xcross[0]; 02096 ycross[1] = ycross[0]; 02097 i2cross[1] = i2cross[0]; 02098 status = 2; 02099 } 02100 xcross[0] = xintersect; 02101 ycross[0] = yintersect; 02102 i2cross[0] = i2; 02103 } 02104 else if ( ncross == 2 && ( count_crossings == 2 || ( 02105 ( ifxsort && ifascend && xintersect < xcross[1] ) || 02106 ( !ifxsort && ifascend && yintersect < ycross[1] ) || 02107 ( ifxsort && !ifascend && xintersect >= xcross[1] ) || 02108 ( !ifxsort && !ifascend && yintersect >= ycross[1] ) ) ) ) 02109 { 02110 xcross[1] = xintersect; 02111 ycross[1] = yintersect; 02112 i2cross[1] = i2; 02113 status = 2; 02114 } 02115 } 02116 } 02117 } 02118 i2m1 = i2; 02119 } 02120 return status; 02121 }