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