From: <airwin@us...>  20090416 17:13:46

Revision: 9807 http://plplot.svn.sourceforge.net/plplot/?rev=9807&view=rev Author: airwin Date: 20090416 17:13:42 +0000 (Thu, 16 Apr 2009) Log Message:  Deal with corner cases in plfill and plfill3 where it was possible to overflow arrays. Thanks to Mark de Wever for discovering the issue. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20090416 14:24:04 UTC (rev 9806) +++ trunk/src/plfill.c 20090416 17:13:42 UTC (rev 9807) @@ 75,7 +75,7 @@ } if (x[0] != x[n1]  y[0] != y[n1]) {  n++; + if(n < PL_MAXPOLY) n++; xpoly[n1] = plP_wcpcx(x[0]); ypoly[n1] = plP_wcpcy(y[0]); } @@ 123,8 +123,8 @@ tx[i] = x[i]; ty[i] = y[i]; tz[i] = z[i]; } if (tx[0] != tx[n1]  ty[0] != ty[n1]  tz[0] != tz[n1]) {  tx[n] = tx[0]; ty[n] = ty[0]; tz[n] = tz[0];  n++; + if(n < PL_MAXPOLY) n++; + tx[n1] = tx[0]; ty[n1] = ty[0]; tz[n1] = tz[0]; } V[0] = tx; V[1] = ty; V[2] = tz; n = plP_clip_poly(n, V, 0, 1, xmin); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <arjenmarkus@us...>  20091208 08:05:16

Revision: 10711 http://plplot.svn.sourceforge.net/plplot/?rev=10711&view=rev Author: arjenmarkus Date: 20091208 08:05:06 +0000 (Tue, 08 Dec 2009) Log Message:  Solve a polygon filling problem that was exhibited in example 25. While it is possible that more such problems exist (missing cases in the filling routine), this commit at least ensures that for not too complicated polygons filling and clipping are done correctly. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091208 06:51:27 UTC (rev 10710) +++ trunk/src/plfill.c 20091208 08:05:06 UTC (rev 10711) @@ 603,7 +603,6 @@ /* Limit case  all vertices are outside of bounding box. So just fill entire * box, *if* the bounding box is completely encircled. */  if ( iclp == 0 ) { if ( inside_lb ) @@ 621,7 +620,6 @@ free( yclp ); }  return; } } @@ 813,8 +811,7 @@ * enough, apparently the crossed_*1 and crossed_*2 variables * are not quite what I expected. */  if ( crossed_left + crossed_right + crossed_down + crossed_up == 1 &&  inside_lb + inside_rb + inside_lu + inside_ru == 4 ) + if ( inside_lb + inside_rb + inside_lu + inside_ru == 4 ) { int dir = circulation( x, y, npts ); PLINT xlim[4], ylim[4]; @@ 825,18 +822,83 @@ xlim[1] = xmax; ylim[1] = ymin; xlim[2] = xmax; ylim[2] = ymax; xlim[3] = xmin; ylim[3] = ymax;  if ( dir > 0 ) + + if ( crossed_left + crossed_right + crossed_down + crossed_up == 1 ) {  incr = 1;  insert = 0 * crossed_left + 1 * crossed_down + 2 * crossed_right +  3 * crossed_up; + if ( dir > 0 ) + { + incr = 1; + insert = 0 * crossed_left + 1 * crossed_down + 2 * crossed_right + + 3 * crossed_up; + } + else + { + incr = 1; + insert = 3 * crossed_left + 2 * crossed_up + 1 * crossed_right + + 0 * crossed_down; + } }  else + + if ( crossed_left + crossed_right == 2 && crossed_down + crossed_up == 0 ) {  incr = 1;  insert = 3 * crossed_left + 2 * crossed_up + 1 * crossed_right +  0 * crossed_down; + if ( xclp[iclp1] == xmin ) + { + if ( dir == 1 ) + { + incr = 1; + insert = 0; + } + else + { + incr = 1; + insert = 3; + } + } + else + { + if ( dir == 1 ) + { + incr = 1; + insert = 1; + } + else + { + incr = 1; + insert = 2; + } + } } + + if ( crossed_left + crossed_right == 0 && crossed_down + crossed_up == 2 ) + { + if ( yclp[iclp1] == ymin ) + { + if ( dir == 1 ) + { + incr = 1; + insert = 1; + } + else + { + incr = 1; + insert = 0; + } + } + else + { + if ( dir == 1 ) + { + incr = 1; + insert = 3; + } + else + { + incr = 1; + insert = 2; + } + } + } + for ( i = 0; i < 4; i++ ) { xclp[iclp] = xlim[insert]; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <andrewross@us...>  20091209 20:35:18

Revision: 10713 http://plplot.svn.sourceforge.net/plplot/?rev=10713&view=rev Author: andrewross Date: 20091209 20:35:05 +0000 (Wed, 09 Dec 2009) Log Message:  Apply style formatting. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091209 08:56:03 UTC (rev 10712) +++ trunk/src/plfill.c 20091209 20:35:05 UTC (rev 10713) @@ 841,7 +841,7 @@ if ( crossed_left + crossed_right == 2 && crossed_down + crossed_up == 0 ) {  if ( xclp[iclp1] == xmin ) + if ( xclp[iclp  1] == xmin ) { if ( dir == 1 ) { @@ 851,7 +851,7 @@ else { incr = 1;  insert = 3; + insert = 3; } } else @@ 864,14 +864,14 @@ else { incr = 1;  insert = 2; + insert = 2; } } } if ( crossed_left + crossed_right == 0 && crossed_down + crossed_up == 2 ) {  if ( yclp[iclp1] == ymin ) + if ( yclp[iclp  1] == ymin ) { if ( dir == 1 ) { @@ 881,7 +881,7 @@ else { incr = 1;  insert = 0; + insert = 0; } } else @@ 894,7 +894,7 @@ else { incr = 1;  insert = 2; + insert = 2; } } } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091216 00:44:44

Revision: 10721 http://plplot.svn.sourceforge.net/plplot/?rev=10721&view=rev Author: airwin Date: 20091216 00:44:29 +0000 (Wed, 16 Dec 2009) Log Message:  Move back to the historical PLINT (as opposed to PLFLT) implementation of pointinpolygon from the current PLINT wrapper for plP_pointinpolygon. The current PLFLT implementation of plP_pointinpolygon is replaced by a PLFLT wrapper for pointinpolygon. I have done these changes because fill_intersection_polygon uses pointinpolygon, and that turned up a bug in the PLFLT form (a test point identical to a polygon vertex returned true). Therefore, I wanted to test and debug using the exact historical PLINT form (rather than a PLINT wrapper for the PLFLT form) since that exact version is the one that has received the most testing in the past. plP_pointinpolygon is currently only used by the software fallback gradient "defined region" approach, and small integer truncation errors would not make much difference in that case. Furthermore, I have plans to replace that approach with one based on fill_intersection_polygon so that in the future it may be possible to simply drop the current plP_pointinpolygon PLFLT wrapper for pointinpolygon completely. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091215 22:28:41 UTC (rev 10720) +++ trunk/src/plfill.c 20091216 00:44:29 UTC (rev 10721) @@ 300,7 +300,7 @@ * Utility functions \**/ static void +void tran( PLINT *a, PLINT *b, PLFLT c, PLFLT d ) { PLINT ta, tb; @@ 312,7 +312,7 @@ *b = (PLINT) floor((double) ( tb * c  ta * d + 0.5 )); } static void +void buildlist( PLINT xp1, PLINT yp1, PLINT xp2, PLINT yp2, PLINT xp3, PLINT yp3, PLINT dinc ) { @@ 363,7 +363,7 @@ } } static void +void addcoord( PLINT xp1, PLINT yp1 ) { PLINT *temp; @@ 385,7 +385,7 @@ buffer[bufferleng++] = yp1; } static int +int compar( const void *pnum1, const void *pnum2 ) { const struct point *pnt1, *pnt2; @@ 991,7 +991,7 @@ * are used to avoid overflow. \**/ static int +int circulation( PLINT *x, PLINT *y, PLINT npts ) { PLFLT xproduct; @@ 1016,37 +1016,33 @@ return direction; } /**\  * int pointinpolygon()  *  * PLINT wrapper for plP_pointinpolygon.  \**/ static int pointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) +/* PLFLT wrapper for pointinpolygon. */ +int +plP_pointinpolygon( PLINT n, const PLFLT *x, const PLFLT *y, PLFLT xp, PLFLT yp ) { int i, return_value;  PLFLT *xflt, *yflt;  if (( xflt = (PLFLT *) malloc( n * sizeof ( PLFLT ))) == NULL ) + PLINT *xint, *yint; + if (( xint = (PLINT *) malloc( n * sizeof ( PLINT ))) == NULL ) {  plexit( "pointinpolygon: Insufficient memory" ); + plexit( "PlP_pointinpolygon: Insufficient memory" ); }  if (( yflt = (PLFLT *) malloc( n * sizeof ( PLFLT ))) == NULL ) + if (( yint = (PLINT *) malloc( n * sizeof ( PLINT ))) == NULL ) {  plexit( "pointinpolygon: Insufficient memory" ); + plexit( "PlP_pointinpolygon: Insufficient memory" ); } for ( i = 0; i < n; i++ ) {  xflt[i] = (PLFLT) x[i];  yflt[i] = (PLFLT) y[i]; + xint[i] = (PLINT) x[i]; + yint[i] = (PLINT) y[i]; }  return_value = plP_pointinpolygon( n, xflt, yflt, (PLFLT) xp, (PLFLT) yp );  free( xflt );  free( yflt ); + return_value = pointinpolygon( n, xint, yint, (PLINT) xp, (PLINT) yp ); + free( xint ); + free( yint ); return return_value; } /**\  * int plP_pointinpolygon() + * int pointinpolygon() * * Returns 1 if the point is inside the polygon, 0 otherwise * Notes: @@ 1058,14 +1054,17 @@ \**/ int plP_pointinpolygon( PLINT n, const PLFLT *x, const PLFLT *y, PLFLT xp, PLFLT yp ) +pointinpolygon( int n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) { int i; int count_crossings;  PLFLT x1, y1, x2, y2, xout, yout, xmax; + PLFLT x1, y1, x2, y2, xpp, ypp, xout, yout, xmax; PLFLT xvp, yvp, xvv, yvv, xv1, yv1, xv2, yv2; PLFLT inprod1, inprod2; + xpp = (PLFLT) xp; + ypp = (PLFLT) yp; + count_crossings = 0; @@ 1090,22 +1089,25 @@ /* Determine for each side whether the line segment between * our two points crosses the vertex */  xvp = xp  xout;  yvp = yp  yout; + xpp = (PLFLT) xp; + ypp = (PLFLT) yp; + xvp = xpp  xout; + yvp = ypp  yout; + for ( i = 0; i < n; i++ ) {  x1 = x[i];  y1 = y[i]; + x1 = (PLFLT) x[i]; + y1 = (PLFLT) y[i]; if ( i < n  1 ) {  x2 = x[i + 1];  y2 = y[i + 1]; + x2 = (PLFLT) x[i + 1]; + y2 = (PLFLT) y[i + 1]; } else {  x2 = x[0];  y2 = y[0]; + x2 = (PLFLT) x[0]; + y2 = (PLFLT) y[0]; } /* Skip zerolength segments */ @@ 1129,11 +1131,11 @@ } /* Line through the two vertices:  * Are xout and xp on either side? */ + * Are xout and xpp on either side? */ xvv = x2  x1; yvv = y2  y1;  xv1 = xp  x1;  yv1 = yp  y1; + xv1 = xpp  x1; + yv1 = ypp  y1; xv2 = xout  x1; yv2 = yout  y1; inprod1 = xv1 * yvv  yv1 * xvv; @@ 1153,7 +1155,6 @@ return ( count_crossings % 2 ); }  /* Fill intersection of two simple (not selfintersecting) polygons. * There must be an even number of edge intersections between the two * polygons (ignoring vertex intersections which touch, but do not cross). This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091218 01:47:46

Revision: 10723 http://plplot.svn.sourceforge.net/plplot/?rev=10723&view=rev Author: airwin Date: 20091218 01:47:34 +0000 (Fri, 18 Dec 2009) Log Message:  Properly scale plP_pointinpolygon PLFLT wrapper for pointinpolygon. Improve notintersect code. Implement a new conditionally compiled version of notpointinpolygon that uses the notintersect code. However, just as a temporary measure, turn off that new version of notpointinpolygon because it gives worse results than the old version. The result is that by default x25c gives good results for the fill and native gradient cases, and also gives pretty good results for the software fallback gradient case. Once I sort out the issue in the notintersect code, I will switch to using the new notinpolygon code. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091216 19:54:02 UTC (rev 10722) +++ trunk/src/plfill.c 20091218 01:47:34 UTC (rev 10723) @@ 2,7 +2,8 @@ * * Polygon pattern fill. *  * Copyright (C) 2004 Alan W. Irwin + * Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 Alan W. Irwin + * Copyright (C) 2005, 2006, 2007, 2008, 2009 Arjen Markus * * This file is part of PLplot. * @@ 50,13 +51,13 @@ buildlist( PLINT, PLINT, PLINT, PLINT, PLINT, PLINT, PLINT ); static int pointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ); +notpointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ); static int circulation( PLINT *x, PLINT *y, PLINT npts ); static void fill_intersection_polygon( PLINT recursion_depth, PLINT ifclip, +fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, void ( *fill )( short *, short *, PLINT ), const PLINT *x1, const PLINT *y1, PLINT i1start, PLINT n1, @@ 64,9 +65,9 @@ const PLINT *if2, PLINT n2 ); static int ifnotintersect( PLINT *xintersect, PLINT *yintersect,  PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2,  PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); +notintersect( PLINT *xintersect, PLINT *yintersect, + PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, + PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); /**\ * void plfill() @@ 489,10 +490,10 @@ plexit( "plP_plfclp: Insufficient memory" ); } }  inside_lb = pointinpolygon( npts, x, y, xmin, ymin );  inside_lu = pointinpolygon( npts, x, y, xmin, ymax );  inside_rb = pointinpolygon( npts, x, y, xmax, ymin );  inside_ru = pointinpolygon( npts, x, y, xmax, ymax ); + inside_lb = !notpointinpolygon( npts, x, y, xmin, ymin ); + inside_lu = !notpointinpolygon( npts, x, y, xmin, ymax ); + inside_rb = !notpointinpolygon( npts, x, y, xmax, ymin ); + inside_ru = !notpointinpolygon( npts, x, y, xmax, ymax ); for ( i = 0; i < npts  1; i++ ) { @@ 1017,12 +1018,13 @@ } /* PLFLT wrapper for pointinpolygon. */ +/* PLFLT wrapper for !notpointinpolygon. */ int plP_pointinpolygon( PLINT n, const PLFLT *x, const PLFLT *y, PLFLT xp, PLFLT yp ) { int i, return_value; PLINT *xint, *yint; + PLFLT xmaximum = fabs( xp ), ymaximum = fabs( yp ), xscale, yscale; if (( xint = (PLINT *) malloc( n * sizeof ( PLINT ))) == NULL ) { plexit( "PlP_pointinpolygon: Insufficient memory" ); @@ 1033,29 +1035,39 @@ } for ( i = 0; i < n; i++ ) {  xint[i] = (PLINT) x[i];  yint[i] = (PLINT) y[i]; + xmaximum = MAX( xmaximum, fabs( x[i] )); + ymaximum = MAX( ymaximum, fabs( y[i] )); }  return_value = pointinpolygon( n, xint, yint, (PLINT) xp, (PLINT) yp ); + xscale = 2.e9 / xmaximum; + yscale = 2.e9 / ymaximum; + for ( i = 0; i < n; i++ ) + { + xint[i] = (PLINT) ( xscale * x[i] ); + yint[i] = (PLINT) ( yscale * y[i] ); + } + return_value = !notpointinpolygon( n, xint, yint, + (PLINT) ( xscale * xp ), (PLINT) ( yscale * yp )); free( xint ); free( yint ); return return_value; } /**\  * int pointinpolygon() + * int notpointinpolygon() *  * Returns 1 if the point is inside the polygon, 0 otherwise + * Returns 0, 1, or 2 depending on whether the test point is definitely + * inside, near the border, or definitely outside the polygon. * Notes:  * Points on the polygon are considered to be outside. * This "Ray casting algorithm" has been described in * http://en.wikipedia.org/wiki/Point_in_polygon. * Logic still needs to be inserted to take care of the "ray passes * through vertex" problem in a numerically robust way. \**/ +#define TEMPORARY_NOTPOINTINPOLYGON_CODE int pointinpolygon( int n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) +notpointinpolygon( int n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) { +#ifdef TEMPORARY_NOTPOINTINPOLYGON_CODE int i; int count_crossings; PLFLT x1, y1, x2, y2, xpp, ypp, xout, yout, xmax; @@ 1153,71 +1165,68 @@ /* Return the result: an even number of crossings means the * point is outside the polygon */  return ( count_crossings % 2 ); + return !( count_crossings % 2 ); } +#else + int i, im1; + int count_crossings = 0; + PLINT xmin, xout, yout, xintersect, yintersect; + + + /* Determine a point outside the polygon */ + + xmin = x[0]; + xout = x[0]; + yout = y[0]; + for ( i = 1; i < n; i++ ) + { + xout = MAX( xout, x[i] ); + xmin = MIN( xmin, x[i] ); + } + /* + 10 to make sure completely outside. */ + xout = xout + ( xout  xmin ) + 10; + + /* Determine whether the line between (xout, yout) and (xp, yp) intersects + * one of the polygon segments. */ + + im1 = n  1; + for ( i = 0; i < n; i++ ) + { + if ( !( x[im1] == x[i] && y[im1] == y[i] ) && + !notintersect( &xintersect, &yintersect, + x[im1], y[im1], x[i], y[i], + xp, yp, xout, yout )) + count_crossings++; + im1 = i; + } + + /* Return the result: an even number of crossings means the + * point is outside the polygon */ + + return !( count_crossings % 2 ); +} +#endif /* TEMPORARY_NOTPOINTINPOLYGON_CODE */ /* Fill intersection of two simple (not selfintersecting) polygons. * There must be an even number of edge intersections between the two * polygons (ignoring vertex intersections which touch, but do not cross). * Eliminate those intersection pairs by recursion (calling the same  * routine again with either the first or second polygon split between  * the two intersecting edges into two independent second polygons.) + * routine twice again with the second polygon split at a boundary defined + * by the first intersection point, all polygon 1 vertices between + * the intersections, and the second intersection point). * Once the recursion has eliminated all intersecting edges, fill or * not using the appropriate polygon depending on whether the first * and second polygons are identical or whether one of them is  * entirely inside the other of them. If ifclip is true, the fill + * entirely inside the other of them. If ifextrapolygon is true, the fill * step will consist of another recursive call to the routine with  * ifclip false and the second polygon set to the clip rectangle. + * ifextrapolygon false, and the second polygon set to an additional + * polygon defined by the stream (not yet implemented). * N.B. it is the calling routine's responsibility to insure the two * polygons are not selfintersecting and do not have duplicate points. */ /* Consider the intersections between the first and second  * polygons (including the possibility of multiple intersections per edge of  * the first polygon). The total number of such intersections must be even.  * Solve for the first pair of intersections (if any) and simplify the problem  * by recursively calling the same routine twice with the same first polygon  * and two modified second polygons consisting of the lefthand segment of the  * second polygon + two intersection points and the righthand segment of the  * second polygon + two intersection points. Note that for the recursive call  * the vertices that were intersection points should not generate new  * intersections with the first unchanging polygon because points on the  * perimeter of polygons are always considered to be outside it. However, due  * to integer truncation and other numerical errors I doubt we can guaranteed  * that numerically in all cases. So we set up a logical array that drops  * segments from the intersection test that are already bounded by  * intersections generated higher in the recursion stack. The result  * of this recursion is deep in the recursion stack you finally end up  * considering the case of the original first polygon and a modified second  * polygon that does not intersect with it at all.  *  *  * . Solve that simple case. Check for identical polygons (taking into account  * the possibilities of different starting points and different directions  * traversing the perimeter), and if so, do a fill of the first polygon (to be  * specific) in that case. If not identical, then do one testinpolygon call to  * decide if the second polygon is totally inside the first. If so, fill the  * second polygon. If not, use one additional testinpolygon call to decide if  * the first polygon is totally inside the second, and if so, fill the first  * polygon.  *  * . Carry a clip flag as part of the recursion, and for the simple case above  * with the flag true, start the recursion all over again (with the clip flag  * false) with the "filling" polygon as the first polygon, and the clip  * rectangle as the the second.  *  *  * Find innermost polygon of (x1,y1) and (x2, y2) for the purposes of  * filling. The number of vertices is returned, and the resulting  * inner (xinner, yinner) polygon calculated for that number of  * vertices. A zero is returned as an error code if the number of  * inner polygon vertices exceeds max_ninner (which should be set  * externally to the size allocated externally for xinner, yinner). A  * zero is also returned if n1 < 3 or (n2 < 3 && n2 != 2). (Note  * that n2 = 2 signals the special case when the second polygon is a  * rectangle.) */ #define MAX_RECURSION_DEPTH 10 void fill_intersection_polygon( PLINT recursion_depth, PLINT ifclip, +fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, void ( *fill )( short *, short *, PLINT ), const PLINT *x1, const PLINT *y1, PLINT i1start, PLINT n1, @@ 1306,9 +1315,9 @@ * *  1 1 * 1 2 1 1 ...  * x + * X * 1  * x + * X * 2 * 1 1 * 1 @@ 1317,11 +1326,11 @@ * 2???2 * *  * "1" marks polygon 1 vertices, "2" marks polygon 2 vertices, "x" + * "1" marks polygon 1 vertices, "2" marks polygon 2 vertices, "X" * marks the intersections, "" stands for part of polygon 1 * that has been previously searched for all possible intersections * from index 0, and "..." means polygon 1 continues  * with more potential intersections both above or below this diagram + * with more potential intersections above and/or below this diagram * before it finally hooks back to connect with the index 0 vertex. * "2???2" stands for parts of polygon 2 that must connect with each other * (since the polygon 1 path between the two intersections is @@ 1359,7 +1368,7 @@ if ( !if2[i2] ) { /* intersect is acted upon only if a series of conditions are met. */  ifintersect = !ifnotintersect( + ifintersect = !notintersect( &xintersect[nintersect], &yintersect[nintersect], x1[i1m1], y1[i1m1], x1[i1], y1[i1], x2[i2m1], y2[i2m1], x2[i2], y2[i2] ); @@ 1500,7 +1509,7 @@ if ( kk >= n2 ) kk = n2; }  fill_intersection_polygon( recursion_depth + 1, ifclip, fill, + fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, x1, y1, i1start, n1, xsplit1, ysplit1, ifsplit1, nsplit1 ); free( xsplit1 ); @@ 1518,7 +1527,7 @@ if ( kk < 0 ) kk += n2; }  fill_intersection_polygon( recursion_depth + 1, ifclip, fill, + fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, x1, y1, i1start, n1, xsplit2, ysplit2, ifsplit2, nsplit2 ); free( xsplit2 ); @@ 1534,20 +1543,28 @@ i1m1 = i1; i1wrap = i1m1; } + if ( nintersect != 0 ) { plwarn( "fill_intersection_polygon: Internal error; nintersect != 0." ); return; }  /* Look for first vertex in Polygon 2 that does not intersect with 1. */ + /* This end phase is reached only if no intersections are found. + * For this case, the intersection of polygon 2 and 1, must be + * either of them (in which case fill with the inner one), or neither + * of them (in which case don't fill at all). */ + + /* Look for first vertex in polygon 2 that is inside polygon 1. */ for ( i2 = 0; i2 < n2; i2++ ) {  if ( !if2[i2] ) + /* Do not bother checking vertices already known to be on the + * boundary with polygon 1. */ + if ( !if2[i2] && !notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] )) break; }  if ( i2 < n2 && pointinpolygon( n1, x1, y1, x2[i2], y2[i2] )) + if ( i2 < n2 ) { /* All of polygon 2 inside polygon 1. */ nfill = n2; @@ 1559,7 +1576,7 @@ /* Look for first vertex in polygon 1 that is inside polygon 2. */ for ( i1 = 0; i1 < n1; i1++ ) {  if ( pointinpolygon( n2, x2, y2, x1[i1], y1[i1] )) + if ( !notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] )) break; } @@ 1571,7 +1588,8 @@ yfiller = y1; } else  plwarn( "fill_intersection_polygon: inscribed polygon not yet implemented" ); + { + } } if ( nfill > 0 ) { @@ 1604,13 +1622,15 @@ * (xintersect, yintersect) via the argument list. */ int ifnotintersect( PLINT * xintersect, PLINT * yintersect,  PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2,  PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ) +notintersect( PLINT * xintersect, PLINT * yintersect, + PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, + PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ) {  PLINT xA2A1, yA2A1, xB2B1, yB2B1;  PLINT xB1A1, yB1A1, xB2A1, yB2A1;  PLFLT factor; + PLFLT factor, fxintersect, fyintersect; + /* These variables are PLFLT not for precision, but to + * avoid integer overflows if they were typed as PLINT. */ + PLFLT xA2A1, yA2A1, xB2B1, yB2B1; + PLFLT xB1A1, yB1A1, xB2A1, yB2A1; /* * Two linear equations to be solved for x and y. * y = ((x  xA1)*yA2 + (xA2  x)*yA1)/(xA2  xA1) @@ 1635,22 +1655,51 @@ yB2B1 = yB2  yB1; factor = xA2A1 * yB2B1  yA2A1 * xB2B1;  /* If two line segments are parallel (including identical) .... */ if ( fabs( factor ) == 0. )  return 1; + { + /* Two line segments are parallel */ + /* For this case check for an intersect for each of the end + * points of A segment against the B segment. */ + fxintersect = xA1; + fyintersect = yA1; + if ( BETW( fxintersect, xB1, xB2 ) && BETW( fyintersect, yB1, yB2 )) + { + *xintersect = fxintersect; + *yintersect = fyintersect; + return 0; + } + else + { + fxintersect = xA2; + fyintersect = yA2; + if ( BETW( fxintersect, xB1, xB2 )  BETW( fyintersect, yB1, yB2 )) + { + *xintersect = fxintersect; + *yintersect = fyintersect; + return 0; + } + else + { + /* parallel and no intersect. */ + return 1; + } + } + } xB1A1 = xB1  xA1; yB1A1 = yB1  yA1; xB2A1 = xB2  xA1; yB2A1 = yB2  yA1; factor = ( xB1A1 * yB2A1  yB1A1 * xB2A1 ) / factor;  *xintersect = factor * xA2A1 + xA1;  *yintersect = factor * yA2A1 + yA1;  /* The x and y range checks (which include end points) are redundant  * with each other for infiniteprecision floatingpoint arithmetic.  * But we don't have that so do these "redundant" checks. */  if ( BETW( *xintersect, xA1, xA2 ) && BETW( *yintersect, yA1, yA2 ) &&  BETW( *xintersect, xB1, xB2 ) && BETW( *yintersect, yB1, yB2 )) + fxintersect = factor * xA2A1 + xA1; + *xintersect = fxintersect; + fyintersect = factor * yA2A1 + yA1; + *yintersect = fyintersect; + + /* The "redundant" x and y segment range checks (which include end points) + * are needed for the vertical and the horizontal cases. */ + if (( BETW( fxintersect, xA1, xA2 )  BETW( fyintersect, yA1, yA2 )) && + ( BETW( fxintersect, xB1, xB2 )  BETW( fyintersect, yB1, yB2 ))) return 0; else return 1; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091218 07:42:53

Revision: 10724 http://plplot.svn.sourceforge.net/plplot/?rev=10724&view=rev Author: airwin Date: 20091218 06:27:40 +0000 (Fri, 18 Dec 2009) Log Message:  Fix scaling for plP_pointinpolygon. (2.e9 should generate integer overflows according to my analysis, but 1.e9 should be okay. However, I had to reduce to 1.e8 for some unknown reason (which I am not going to pursue since plP_pointinpolygon will probably be removed soon in any case). Fix bugs in notintersect. The result of these changes is the new notpointinpolygon code works (as well as the old code) for example 25 so use "#define NEW_NOTPOINTINPOLYGON_CODE" to compile that new code while avoiding compiling the old code (which I will retain for a while longer for reference). Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091218 01:47:34 UTC (rev 10723) +++ trunk/src/plfill.c 20091218 06:27:40 UTC (rev 10724) @@ 1038,8 +1038,8 @@ xmaximum = MAX( xmaximum, fabs( x[i] )); ymaximum = MAX( ymaximum, fabs( y[i] )); }  xscale = 2.e9 / xmaximum;  yscale = 2.e9 / ymaximum; + xscale = 1.e8 / xmaximum; + yscale = 1.e8 / ymaximum; for ( i = 0; i < n; i++ ) { xint[i] = (PLINT) ( xscale * x[i] ); @@ 1063,11 +1063,49 @@ * through vertex" problem in a numerically robust way. \**/ #define TEMPORARY_NOTPOINTINPOLYGON_CODE +#define NEW_NOTPOINTINPOLYGON_CODE int notpointinpolygon( int n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) { #ifdef TEMPORARY_NOTPOINTINPOLYGON_CODE +#ifdef NEW_NOTPOINTINPOLYGON_CODE + int i, im1; + int count_crossings = 0; + PLINT xmin, xout, yout, xintersect, yintersect; + + + /* Determine a point outside the polygon */ + + xmin = x[0]; + xout = x[0]; + yout = y[0]; + for ( i = 1; i < n; i++ ) + { + xout = MAX( xout, x[i] ); + xmin = MIN( xmin, x[i] ); + } + /* + 10 to make sure completely outside. */ + xout = xout + ( xout  xmin ) + 10; + + /* Determine whether the line between (xout, yout) and (xp, yp) intersects + * one of the polygon segments. */ + + im1 = n  1; + for ( i = 0; i < n; i++ ) + { + if ( !( x[im1] == x[i] && y[im1] == y[i] ) && + !notintersect( &xintersect, &yintersect, + x[im1], y[im1], x[i], y[i], + xp, yp, xout, yout )) + count_crossings++; + im1 = i; + } + + /* Return the result: an even number of crossings means the + * point is outside the polygon */ + + return !( count_crossings % 2 ); +} +#else /* NEW_NOTPOINTINPOLYGON_CODE */ int i; int count_crossings; PLFLT x1, y1, x2, y2, xpp, ypp, xout, yout, xmax; @@ 1167,45 +1205,7 @@ return !( count_crossings % 2 ); } #else  int i, im1;  int count_crossings = 0;  PLINT xmin, xout, yout, xintersect, yintersect;    /* Determine a point outside the polygon */   xmin = x[0];  xout = x[0];  yout = y[0];  for ( i = 1; i < n; i++ )  {  xout = MAX( xout, x[i] );  xmin = MIN( xmin, x[i] );  }  /* + 10 to make sure completely outside. */  xout = xout + ( xout  xmin ) + 10;   /* Determine whether the line between (xout, yout) and (xp, yp) intersects  * one of the polygon segments. */   im1 = n  1;  for ( i = 0; i < n; i++ )  {  if ( !( x[im1] == x[i] && y[im1] == y[i] ) &&  !notintersect( &xintersect, &yintersect,  x[im1], y[im1], x[i], y[i],  xp, yp, xout, yout ))  count_crossings++;  im1 = i;  }   /* Return the result: an even number of crossings means the  * point is outside the polygon */   return !( count_crossings % 2 ); } #endif /* TEMPORARY_NOTPOINTINPOLYGON_CODE */ +#endif /* NEW_NOTPOINTINPOLYGON_CODE */ /* Fill intersection of two simple (not selfintersecting) polygons. * There must be an even number of edge intersections between the two * polygons (ignoring vertex intersections which touch, but do not cross). @@ 1672,7 +1672,7 @@ { fxintersect = xA2; fyintersect = yA2;  if ( BETW( fxintersect, xB1, xB2 )  BETW( fyintersect, yB1, yB2 )) + if ( BETW( fxintersect, xB1, xB2 ) && BETW( fyintersect, yB1, yB2 )) { *xintersect = fxintersect; *yintersect = fyintersect; @@ 1698,8 +1698,8 @@ /* The "redundant" x and y segment range checks (which include end points) * are needed for the vertical and the horizontal cases. */  if (( BETW( fxintersect, xA1, xA2 )  BETW( fyintersect, yA1, yA2 )) &&  ( BETW( fxintersect, xB1, xB2 )  BETW( fyintersect, yB1, yB2 ))) + if (( BETW( fxintersect, xA1, xA2 ) && BETW( fyintersect, yA1, yA2 )) && + ( BETW( fxintersect, xB1, xB2 ) && BETW( fyintersect, yB1, yB2 ))) return 0; else return 1; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091219 23:00:24

Revision: 10725 http://plplot.svn.sourceforge.net/plplot/?rev=10725&view=rev Author: airwin Date: 20091219 23:00:16 +0000 (Sat, 19 Dec 2009) Log Message:  Improved handling of cases where the intersection between two line segments is near the ends of the line segments being tested. This change affects the USE_FILL_INTERSECTION_POLYGON ON case (via a call of notcrossed in both fill_intersection_polygon and pointinpolygon) and USE_FILL_INTERSECTION_POLYGON OFF cases (via the call of notcrossed from pointinpolygon). Improved handling of the end phase of fill_intersection_polygon where polygon 1 and the split of polygon 2 do not cross each other. Bug fixes for fill_intersection_polygon. These changes only affect the case where USE_FILL_INTERSECTION_POLYGON is ON. example 25 results continue to look good for the default USE_FILL_INTERSECTION_POLYGON OFF case. However, there are still some obvious issues (not memory management issues according to valgrind) for the experimental USE_FILL_INTERSECTION_POLYGON ON case. For this case, the example 25 results sometimes fill in the wrong locations, and I am currently working on tracking down and fixing all of these issues. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091218 06:27:40 UTC (rev 10724) +++ trunk/src/plfill.c 20091219 23:00:16 UTC (rev 10725) @@ 27,9 +27,42 @@ #define INSIDE( ix, iy ) ( BETW( ix, xmin, xmax ) && BETW( iy, ymin, ymax )) #define DTOR ( PI / 180. ) #define BINC 50 +#define DTOR ( PI / 180. ) +#define BINC 50 +/* Nearborder comparison criterion (NBCC). */ +#define PL_NBCC 2 +/* Variant of BETW that returns true if between or within PL_NBCC of it. */ +#define BETW_NBCC( ix, ia, ib ) ((( ix ) <= ( ia + PL_NBCC ) && ( ix ) >= ( ib  PL_NBCC ))  (( ix ) >= ( ia  PL_NBCC ) && ( ix ) <= ( ib + PL_NBCC ))) +/* Status codes ORed together in the return from notcrossed. + * PL_NOT_CROSSED is set whenever the two line segments definitely + * (i.e., not near the ends or completely apart) do not cross each + * other. + * + * PL_NEAR_A1 is set if the intersect is near (+/ PL_NBCC) the first + * end of line segment A. + * + * PL_NEAR_A2 is set if the intersect is near (+/ PL_NBCC) the second + * end of line segment A. + * + * PL_NEAR_B1 is set if the intersect is near (+/ PL_NBCC) the first + * end of line segment B. + * + * PL_NEAR_B2 is set if the intersect is near (+/ PL_NBCC) the second + * end of line segment B. + * + * PL_PARALLEL is set if the two line segments are parallel with each other. + */ +enum PL_CrossedStatus +{ + PL_NOT_CROSSED = 0x1, + PL_NEAR_A1 = 0x2, + PL_NEAR_A2 = 0x4, + PL_NEAR_B1 = 0x8, + PL_NEAR_B2 = 0x10, + PL_PARALLEL = 0x20 +}; + struct point { PLINT x, y; @@ 65,9 +98,9 @@ const PLINT *if2, PLINT n2 ); static int notintersect( PLINT *xintersect, PLINT *yintersect,  PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2,  PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); +notcrossed( PLINT *xintersect, PLINT *yintersect, + PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, + PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); /**\ * void plfill() @@ 1063,12 +1096,13 @@ * through vertex" problem in a numerically robust way. \**/ +/* Temporary until get rid of old code altogether. */ #define NEW_NOTPOINTINPOLYGON_CODE int notpointinpolygon( int n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp ) { #ifdef NEW_NOTPOINTINPOLYGON_CODE  int i, im1; + int i, im1, ifnotcrossed; int count_crossings = 0; PLINT xmin, xout, yout, xintersect, yintersect; @@ 1092,18 +1126,28 @@ im1 = n  1; for ( i = 0; i < n; i++ ) {  if ( !( x[im1] == x[i] && y[im1] == y[i] ) &&  !notintersect( &xintersect, &yintersect,  x[im1], y[im1], x[i], y[i],  xp, yp, xout, yout ))  count_crossings++; + if ( !( x[im1] == x[i] && y[im1] == y[i] )) + { + ifnotcrossed = notcrossed( &xintersect, &yintersect, + x[im1], y[im1], x[i], y[i], + xp, yp, xout, yout ); + + if ( !ifnotcrossed ) + count_crossings++; + else if ( ifnotcrossed & ( PL_NEAR_A1  PL_NEAR_A2  PL_NEAR_B1  PL_NEAR_B2 )) + return 1; + } im1 = i; }  /* Return the result: an even number of crossings means the  * point is outside the polygon */   return !( count_crossings % 2 ); + /* return 0 if the test point is definitely inside + * (count_crossings odd), return 1 if the test point is near (see + * above logic), and return 2 if the test point is definitely + * outside the border (count_crossings even). */ + if (( count_crossings % 2 ) == 1 ) + return 0; + else + return 2; } #else /* NEW_NOTPOINTINPOLYGON_CODE */ int i; @@ 1236,12 +1280,13 @@ PLINT i1, i1m1, i1wrap, i1wraplast, i2, i2m1, i2wrap, i2wraplast, kk, kkstart1, kkstart21, kkstart22,  k, kstart, range1, range21, range22, nintersect, + k, kstart, range1, range21, range22, ncrossed, nsplit1, nsplit2;  PLINT xintersect[2], yintersect[2], ifintersect; + PLINT xintersect[2], yintersect[2], ifcrossed; PLINT *xsplit1, *ysplit1, *ifsplit1, *xsplit2, *ysplit2, *ifsplit2;  PLINT ifill, nfill = 0; + PLINT ifill, nfill = 0, + ifnotpolygon1inpolygon2, ifnotpolygon2inpolygon1; const PLINT *xfiller, *yfiller; short *xfill, *yfill; @@ 1271,7 +1316,7 @@ } i1m1 = i1start  1;  if ( i1m1 <= 0 ) + if ( i1m1 < 0 ) i1m1 = n1  1; for ( i1 = i1start; i1 < n1; i1++ ) @@ 1297,7 +1342,7 @@ if ( i2 < n2 ) {  plwarn( "fill_intersection_polygon: Internal error; i < n2." ); + plwarn( "fill_intersection_polygon: Internal error; i2 < n2." ); return; } @@ 1356,9 +1401,12 @@ * intersections to be found between polygon 1 and some polygon 2 * split, and in that case we move on to the end phase below. */  nintersect = 0;  i1m1 = n1  1;  i1wrap = 1; + ncrossed = 0; + i1wrap = i1start  1; + if ( i1wrap < 0 ) + i1m1 = n1  1; + else + i1m1 = i1wrap; for ( i1 = i1start; i1 < n1; i1++ ) { i2m1 = n2  1; @@ 1367,34 +1415,23 @@ { if ( !if2[i2] ) {  /* intersect is acted upon only if a series of conditions are met. */  ifintersect = !notintersect(  &xintersect[nintersect], &yintersect[nintersect], + /* ifcrossed true only if there is a definite crossing of + * the two line segments with the intersect not being + * near (+/ PL_NBCC) the ends. ToDo. This test + * must be made more elaborate if the polygons definitely + * cross near one or both of their vertices. */ + ifcrossed = !notcrossed( + &xintersect[ncrossed], &yintersect[ncrossed], x1[i1m1], y1[i1m1], x1[i1], y1[i1], x2[i2m1], y2[i2m1], x2[i2], y2[i2] );  ifintersect = ifintersect && !(  xintersect[nintersect] == x1[i1m1] &&  yintersect[nintersect] == y1[i1m1] );  ifintersect = ifintersect && !(  xintersect[nintersect] == x1[i1] &&  yintersect[nintersect] == y1[i1] );  ifintersect = ifintersect && !(  xintersect[nintersect] == x2[i2m1] &&  yintersect[nintersect] == y2[i2m1] );  ifintersect = ifintersect && !(  xintersect[nintersect] == x2[i2] &&  yintersect[nintersect] == y2[i2] );  if ( ifintersect ) + if ( ifcrossed ) {  /* The above test must be more elaborate if the intersection  * point is a vertex of the first polygon, the second polygon  * or both. ToDo! */  if ( nintersect == 0 ) + if ( ncrossed == 0 ) { i1wraplast = i1wrap; i2wraplast = i2wrap;  nintersect++; + ncrossed++; } else { @@ 1544,53 +1581,81 @@ i1wrap = i1m1; }  if ( nintersect != 0 ) + if ( ncrossed != 0 ) {  plwarn( "fill_intersection_polygon: Internal error; nintersect != 0." ); + plwarn( "fill_intersection_polygon: Internal error; ncrossed != 0." ); return; } /* This end phase is reached only if no intersections are found.  * For this case, the intersection of polygon 2 and 1, must be + * For this case, there are 5 valid and 1 invalid combination of + * the following conditions:the intersection of polygon 2 and 1, must be * either of them (in which case fill with the inner one), or neither * of them (in which case don't fill at all). */  /* Look for first vertex in polygon 2 that is inside polygon 1. */ + /* Classify polygon 1 by looking for first vertex in polygon 1 + * that is definitely inside or outside polygon 2. */ + for ( i1 = 0; i1 < n1; i1++ ) + { + if (( ifnotpolygon1inpolygon2 = + notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] )) != 1 ) + break; + } + + /* Classify polygon 2 by looking for first vertex in polygon 2 + * that is definitely inside or outside polygon 1. */ + ifnotpolygon2inpolygon1 = 1; for ( i2 = 0; i2 < n2; i2++ ) { /* Do not bother checking vertices already known to be on the * boundary with polygon 1. */  if ( !if2[i2] && !notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] )) + if ( !if2[i2] && ( ifnotpolygon2inpolygon1 = + notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] )) != 1 ) break; }  if ( i2 < n2 ) + if ( ifnotpolygon2inpolygon1 == 0 ) {  /* All of polygon 2 inside polygon 1. */  nfill = n2;  xfiller = x2;  yfiller = y2; + /* Polygon 2 definitely inside polygon 1. */ + if ( ifnotpolygon1inpolygon2 != 0 ) + { + nfill = n2; + xfiller = x2; + yfiller = y2; + } + else + plwarn( "fill_intersection_polygon: Internal error; no intersections found but each polygon definitely inside the other!" ); }  else + else if ( ifnotpolygon2inpolygon1 == 1 ) {  /* Look for first vertex in polygon 1 that is inside polygon 2. */  for ( i1 = 0; i1 < n1; i1++ ) + /* Polygon 2 vertices near polygon 1 border. */ + if ( ifnotpolygon1inpolygon2 != 0 ) {  if ( !notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] ))  break; + /* Polygon 1 vertices near polygon 2 border or definitely outside it. */ + nfill = n2; + xfiller = x2; + yfiller = y2; }   if ( i1 < n1 ) + else {  /* All of polygon 1 inside polygon 2 */ + /* Polygon 1 vertices definitely inside polygon 2 border. */ nfill = n1; xfiller = x1; yfiller = y1; }  else + } + else + { + /* Polygon 2 vertices definitely outside polygon 1 border. */ + if ( ifnotpolygon1inpolygon2 != 2 ) { + /* Polygon 1 vertices definitely inside polygon 2 border. */ + nfill = n1; + xfiller = x1; + yfiller = y1; } } + if ( nfill > 0 ) { if (( xfill = (short *) malloc( nfill * sizeof ( short ))) == NULL ) @@ 1614,23 +1679,27 @@ return; } /* Find if an intersection (N.B. including end points) exists between  * straight line segments A and B defined by their endpoints xA1, yA1,  * xA2, yA2, xB1, yB1, xB2, yB2. If not return true (1 if parallel,  * 1 if intersecting away from segments). If true, return false.  * Except for the parallel line case always return the intersection,  * (xintersect, yintersect) via the argument list. */ +/* Returns a 0 status code if the two line segments A, and B defined + * by their end points (xA1, yA1, xA2, yA2, xB1, yB1, xB2, and yB2) + * definitely (i.e., intersection point is not near the ends of either + * of the line segments) cross each other. Otherwise, return a status + * code specifying the type of noncrossing (i.e., completely + * separate, near one of the ends, parallel). Return the actual + * intersection (or average of the x end points and y end points for + * the parallel, not crossed case) via the argument list pointers to + * xintersect and yintersect. */ int notintersect( PLINT * xintersect, PLINT * yintersect,  PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2,  PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ) +notcrossed( PLINT * xintersect, PLINT * yintersect, + PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, + PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ) { PLFLT factor, fxintersect, fyintersect; /* These variables are PLFLT not for precision, but to * avoid integer overflows if they were typed as PLINT. */ PLFLT xA2A1, yA2A1, xB2B1, yB2B1; PLFLT xB1A1, yB1A1, xB2A1, yB2A1; + PLINT status = 0; /* * Two linear equations to be solved for x and y. * y = ((x  xA1)*yA2 + (xA2  x)*yA1)/(xA2  xA1) @@ 1658,49 +1727,60 @@ if ( fabs( factor ) == 0. ) { /* Two line segments are parallel */  /* For this case check for an intersect for each of the end  * points of A segment against the B segment. */  fxintersect = xA1;  fyintersect = yA1;  if ( BETW( fxintersect, xB1, xB2 ) && BETW( fyintersect, yB1, yB2 )) + status = status  PL_PARALLEL; + /* Choice of intersect is arbitrary in this case. Choose A1 if + * that lies near or in B. Otherwise, choose A2 if that lies near + * or in B. Otherwise, choose the average point. */ + if (( BETW_NBCC( xA1, xB1, xB2 ) && BETW_NBCC( yA1, yB1, yB2 ))) {  *xintersect = fxintersect;  *yintersect = fyintersect;  return 0; + fxintersect = xA1; + fyintersect = yA1; }  else + else if (( BETW_NBCC( xA2, xB1, xB2 ) && BETW_NBCC( yA2, yB1, yB2 ))) { fxintersect = xA2; fyintersect = yA2;  if ( BETW( fxintersect, xB1, xB2 ) && BETW( fyintersect, yB1, yB2 ))  {  *xintersect = fxintersect;  *yintersect = fyintersect;  return 0;  }  else  {  /* parallel and no intersect. */  return 1;  } } + else + { + fxintersect = 0.25 * ( xA1 + xA2 + xB1 + xB2 ); + fyintersect = 0.25 * ( yA1 + yA2 + yB1 + yB2 ); + } }  xB1A1 = xB1  xA1;  yB1A1 = yB1  yA1;  xB2A1 = xB2  xA1;  yB2A1 = yB2  yA1; + else + { + xB1A1 = xB1  xA1; + yB1A1 = yB1  yA1; + xB2A1 = xB2  xA1; + yB2A1 = yB2  yA1;  factor = ( xB1A1 * yB2A1  yB1A1 * xB2A1 ) / factor;  fxintersect = factor * xA2A1 + xA1; + factor = ( xB1A1 * yB2A1  yB1A1 * xB2A1 ) / factor; + fxintersect = factor * xA2A1 + xA1; + fyintersect = factor * yA2A1 + yA1; + } + /* The "redundant" x and y segment range checks (which include near the + * end points) are needed for the vertical and the horizontal cases. */ + if (( BETW_NBCC( fxintersect, xA1, xA2 ) && BETW_NBCC( fyintersect, yA1, yA2 )) && + ( BETW_NBCC( fxintersect, xB1, xB2 ) && BETW_NBCC( fyintersect, yB1, yB2 ))) + { + /* The intersect is close (within +/ PL_NBCC) to an end point or + * corresponds to a definite crossing of the two line segments. + * Find out which. */ + if ( fabs( fxintersect  xA1 ) <= PL_NBCC && fabs( fyintersect  yA1 ) <= PL_NBCC ) + status = status  PL_NEAR_A1; + else if ( fabs( fxintersect  xA2 ) <= PL_NBCC && fabs( fyintersect  yA2 ) <= PL_NBCC ) + status = status  PL_NEAR_A2; + else if ( fabs( fxintersect  xB1 ) <= PL_NBCC && fabs( fyintersect  yB1 ) <= PL_NBCC ) + status = status  PL_NEAR_B1; + else if ( fabs( fxintersect  xB2 ) <= PL_NBCC && fabs( fyintersect  yB2 ) <= PL_NBCC ) + status = status  PL_NEAR_B2; + /* N.B. if none of the above conditions hold then status remains at + * zero to signal we have a definite crossing. */ + } + else + status = status  PL_NOT_CROSSED; *xintersect = fxintersect;  fyintersect = factor * yA2A1 + yA1; *yintersect = fyintersect;  /* The "redundant" x and y segment range checks (which include end points)  * are needed for the vertical and the horizontal cases. */  if (( BETW( fxintersect, xA1, xA2 ) && BETW( fyintersect, yA1, yA2 )) &&  ( BETW( fxintersect, xB1, xB2 ) && BETW( fyintersect, yB1, yB2 )))  return 0;  else  return 1; + return status; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091221 07:39:37

Revision: 10726 http://plplot.svn.sourceforge.net/plplot/?rev=10726&view=rev Author: airwin Date: 20091221 07:39:21 +0000 (Mon, 21 Dec 2009) Log Message:  Deal with the case where polygon 1 _must_ be split rather than polygon 2 in fill_intersection_polygon. This change only affects DUSE_FILL_INTERSECTION_POLYGON=ON results and should make no difference for DUSE_FILL_INTERSECTION_POLYGON=OFF results (the current default). The DUSE_FILL_INTERSECTION_POLYGON=ON result is that page 1 of example 25 is correct other than a missing fill issue for subpage 8. There are additional fill and gradient issues for further pages in example 25 so more debugging work is needed, but this result is the best DUSE_FILL_INTERSECTION_POLYGON=ON result obtained to date. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091219 23:00:16 UTC (rev 10725) +++ trunk/src/plfill.c 20091221 07:39:21 UTC (rev 10726) @@ 1437,6 +1437,79 @@ { /* Have discovered the first two intersections for * polygon 1 at i1 = i1start or above. */ + + /* Calculate polygon 2 index range in split 1 + * (the split that normally proceeds beyond + * the second intersect with ascending i2 + * values). */ + range21 = i2wraplast  i2wrap; + if ( range21 < 0 ) + range21 += n2; + /* Wrap values range between 1 and n2  2 so + * the the smallest untransformed range21 + * value is n2 + 1, the smallest transformed + * range21 value is 0 (only obtained if + * i2wraplast == i2wrap), and the largest + * transformed range21 value is n2  1. */ + if ( range21 == 0 ) + { + /* For this special case the above ascii art + * does not apply, and we must use the + * following diagram instead: + * + *  1 1 ... + * + * + * + * 2???2 X X 2???2 + * + * 1 1 + * + * + * N.B. no valid split of polygon 2 is + * possible with this diagram. However, + * this diagrem can be classified the same + * as the first diagram, but with the roles + * of polygon 1 and polygon 2 swapped so + * that polygon 1 is the one that gets split + * instead of polygon 2. In fact, swapping + * those two polygons (and starting the new + * polygon 1 just at i2 of the old polygon 2 + * to insure the present two intersections + * are immediately (at least for the first + * new polygon segment 1) found and dealt + * with to eliminate the chance of an + * infinite "swap" loop) is the only clean + * way to deal with this situation. */ + if (( xsplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ysplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ifsplit1 = (PLINT *) calloc( n1, sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + kk = i2; + for ( i2 = 0; i2 < n2; i2++ ) + { + xsplit2[kk] = x2[i2]; + ysplit2[kk++] = y2[i2]; + if ( kk == n2 ) + kk = 0; + } + fill_intersection_polygon( + recursion_depth + 1, ifextrapolygon, fill, + xsplit2, ysplit2, 0, n2, + x1, y1, ifsplit1, n1 ); + free( xsplit2 ); + free( ysplit2 ); + free( ifsplit1 ); + return; + } /* New i1start is the largest nonnegative polygon 1 * index below the last detected intersect. */ i1start = MAX( i1wrap, 0 ); @@ 1455,18 +1528,6 @@ * kkstart1 of polygon 1, and the second intersect. */ kkstart1 = i1wraplast + 1;  /* Split 1 of polygon2 consists of the  * boundary + range21 points between kkstart21  * (= i2) and i2wraplast in ascending order of  * polygon 2 indices. N.B. if range21 is zero  * below we change that to n2, i.e., use the  * whole range of polygon 2 in ascending  * order. For this case, range22 (below) will  * be zero and range1 (above) must be nonzero  * (in order to have two intersections). */  range21 = i2wraplast  i2wrap;  if ( range21 <= 0 )  range21 += n2; kkstart21 = i2; /* Split 2 of polygon 2 consists of the @@ 1546,7 +1607,9 @@ if ( kk >= n2 ) kk = n2; }  fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, + + fill_intersection_polygon( + recursion_depth + 1, ifextrapolygon, fill, x1, y1, i1start, n1, xsplit1, ysplit1, ifsplit1, nsplit1 ); free( xsplit1 ); @@ 1564,7 +1627,9 @@ if ( kk < 0 ) kk += n2; }  fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, + + fill_intersection_polygon( + recursion_depth + 1, ifextrapolygon, fill, x1, y1, i1start, n1, xsplit2, ysplit2, ifsplit2, nsplit2 ); free( xsplit2 ); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091222 02:31:12

Revision: 10727 http://plplot.svn.sourceforge.net/plplot/?rev=10727&view=rev Author: airwin Date: 20091222 02:31:06 +0000 (Tue, 22 Dec 2009) Log Message:  Fix ifi2 condition for looking for crossings between polygon 1 and 2. For split2 store data in order of increasing index of polygon 2 to preserve the orientation. (These changes to fill_intersection_polygon made page 2 of example 25, where both polygons have positive orientation (see http://en.wikipedia.org/wiki/Curve_orientation for a discussion of this) give correct filling results for the first time for the DFILL_INTERSECTION_POLYGON=ON case. Transform polygon 1 to have positive orientation in the DFILL_INTERSECTION_POLYGON=ON part of plP_plfclp. This made page 1 of example 25 (where polygon 1 had negative orientation) work correctly. Thanks to Arjen for this idea. There are still issues with pages 3 and above of example 25 for DFILL_INTERSECTION_POLYGON=ON. These changes should not affect the default DFILL_INTERSECTION_POLYGON=OFF case. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091221 07:39:21 UTC (rev 10726) +++ trunk/src/plfill.c 20091222 02:31:06 UTC (rev 10727) @@ 102,6 +102,9 @@ PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ); +static int +positive_orientation( PLINT n, const PLINT *x, const PLINT *y ); + /**\ * void plfill() * @@ 457,22 +460,22 @@ if ( npts < 3  !draw ) return; #ifdef USE_FILL_INTERSECTION_POLYGON  PLINT *x1, *y1, i1start = 0, i, im1, n1; + PLINT *x10, *y10, *x1, *y1, i1start = 0, i, im1, n1, n1m1; PLINT x2[4] = { xmin, xmax, xmax, xmin }; PLINT y2[4] = { ymin, ymin, ymax, ymax }; PLINT if2[4] = { 0, 0, 0, 0 }; PLINT n2 = 4;  if (( x1 = (PLINT *) malloc( npts * sizeof ( PLINT ))) == NULL ) + if (( x10 = (PLINT *) malloc( npts * sizeof ( PLINT ))) == NULL ) { plexit( "plP_plfclp: Insufficient memory" ); }  if (( y1 = (PLINT *) malloc( npts * sizeof ( PLINT ))) == NULL ) + if (( y10 = (PLINT *) malloc( npts * sizeof ( PLINT ))) == NULL ) { plexit( "plP_plfclp: Insufficient memory" ); }  /* Polygon 2 obviously has no dups, but get rid of them in polygon  * 1 if they exist. ToDo: Deal with selfintersecting polygon 1  * case as well. */ + /* Polygon 2 obviously has no dups nor two consective segments that + * are parallel, but get rid of those type of segments in polygon 1 + * if they exist. */ im1 = npts  1; n1 = 0; @@ 480,12 +483,49 @@ { if ( !( x[i] == x[im1] && y[i] == y[im1] )) {  x1[n1] = x[i];  y1[n1++] = y[i]; + x10[n1] = x[i]; + y10[n1++] = y[i]; } im1 = i; } + /* Must have at least three points that satisfy the above criteria. */ + if ( n1 < 3 ) + { + free( x10 ); + free( y10 ); + return; + } + + /* Polygon 2 obviously has a positive orientation (i.e., as you + * ascend in index along the boundary, the points just adjacent to + * the boundary and on the left are interior points for the + * polygon), but enforce this condition demanded by + * fill_intersection_polygon for polygon 1 as well. */ + if ( positive_orientation( n1, x10, y10 )) + { + x1 = x10; + y1 = y10; + } + else + { + if (( x1 = (PLINT *) malloc( n1 * sizeof ( PLINT ))) == NULL ) + { + plexit( "plP_plfclp: Insufficient memory" ); + } + if (( y1 = (PLINT *) malloc( n1 * sizeof ( PLINT ))) == NULL ) + { + plexit( "plP_plfclp: Insufficient memory" ); + } + n1m1 = n1  1; + for ( i = 0; i < n1; i++ ) + { + x1[n1m1  i] = x10[i]; + y1[n1m1  i] = y10[i]; + } + free( x10 ); + free( y10 ); + } fill_intersection_polygon( 0, 0, draw, x1, y1, i1start, n1, x2, y2, if2, n2 ); free( x1 ); free( y1 ); @@ 1250,25 +1290,34 @@ return !( count_crossings % 2 ); } #endif /* NEW_NOTPOINTINPOLYGON_CODE */ /* Fill intersection of two simple (not selfintersecting) polygons.  * There must be an even number of edge intersections between the two  * polygons (ignoring vertex intersections which touch, but do not cross).  * Eliminate those intersection pairs by recursion (calling the same  * routine twice again with the second polygon split at a boundary defined  * by the first intersection point, all polygon 1 vertices between  * the intersections, and the second intersection point).  * Once the recursion has eliminated all intersecting edges, fill or  * not using the appropriate polygon depending on whether the first  * and second polygons are identical or whether one of them is  * entirely inside the other of them. If ifextrapolygon is true, the fill  * step will consist of another recursive call to the routine with + +#define MAX_RECURSION_DEPTH 10 + +/* Fill intersection of two simple (not selfintersecting) polygons + * that both have a positive orientation (see + * http://en.wikipedia.org/wiki/Curve_orientation). That is, as you + * traverse the boundary in index order, the inside area of the + * polygon is always on the left. This requirement simplifies the + * logic of fill_instersection_polygon. N.B. it is the calling + * routine's responsibility to insure the two polygons do not have + * duplicate points, are not selfintersecting, and have positive + * orientation. + * + * Two polygons that do not self intersect must have an even number of + * edge crossings between them. (ignoring vertex intersections which + * touch, but do not cross). fill_intersection_polygon eliminates + * those intersection crossings by recursion (calling the same routine + * twice again with the second polygon split at a boundary defined by + * the first intersection point, all polygon 1 vertices between the + * intersections, and the second intersection point). Once the + * recursion has eliminated all crossing edges, fill or not using the + * appropriate polygon depending on whether the first and second + * polygons are identical or whether one of them is entirely inside + * the other of them. If ifextrapolygon is true, the fill step will + * consist of another recursive call to the routine with * ifextrapolygon false, and the second polygon set to an additional  * polygon defined by the stream (not yet implemented).  * N.B. it is the calling routine's responsibility to insure the two  * polygons are not selfintersecting and do not have duplicate points. */ + * polygon defined by the stream (not yet implemented). */  #define MAX_RECURSION_DEPTH 10 void fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, void ( *fill )( short *, short *, PLINT ), @@ 1281,8 +1330,8 @@ i2, i2m1, i2wrap, i2wraplast, kk, kkstart1, kkstart21, kkstart22, k, kstart, range1, range21, range22, ncrossed,  nsplit1, nsplit2;  PLINT xintersect[2], yintersect[2], ifcrossed; + nsplit1, nsplit2, nsplit2m1; + PLINT xintersect[2], yintersect[2], ifnotcrossed; PLINT *xsplit1, *ysplit1, *ifsplit1, *xsplit2, *ysplit2, *ifsplit2; PLINT ifill, nfill = 0, @@ 1413,18 +1462,14 @@ i2wrap = 1; for ( i2 = 0; i2 < n2; i2++ ) {  if ( !if2[i2] ) + if ( !( if2[i2] && if2[i2m1] )) {  /* ifcrossed true only if there is a definite crossing of  * the two line segments with the intersect not being  * near (+/ PL_NBCC) the ends. ToDo. This test  * must be made more elaborate if the polygons definitely  * cross near one or both of their vertices. */  ifcrossed = !notcrossed( + ifnotcrossed = notcrossed( &xintersect[ncrossed], &yintersect[ncrossed], x1[i1m1], y1[i1m1], x1[i1], y1[i1], x2[i2m1], y2[i2m1], x2[i2], y2[i2] );  if ( ifcrossed ) + /* Use only definite crossing case. */ + if ( !ifnotcrossed ) { if ( ncrossed == 0 ) { @@ 1501,6 +1546,12 @@ if ( kk == n2 ) kk = 0; } + /* N.B. the positive orientation of split2 + * is preserved since the index order is + * the same as that of polygon 2, and by + * assumption that polygon and polygon 1 + * have identical positive + * orientations. */ fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, xsplit2, ysplit2, 0, n2, @@ 1567,33 +1618,41 @@ plexit( "fill_intersection_polygon: Insufficient memory" ); } /* Common boundary between split1 and split2. */  k = 0;  xsplit1[k] = xintersect[0];  ysplit1[k] = yintersect[0];  ifsplit1[k] = 1;  xsplit2[k] = xintersect[0];  ysplit2[k] = yintersect[0];  ifsplit2[k] = 1;  kstart = k + 1;  kk = kkstart1; + /* N.B. Although basic index arithmetic for + * split 2 is done in negative orientation + * order because the index is decrementing + * relative to the index of split 2, actually + * store results in reverse order to preserve + * the positive orientation that by assumption + * both polygon 1 and 2 have. */ + k = 0; + xsplit1[k] = xintersect[0]; + ysplit1[k] = yintersect[0]; + ifsplit1[k] = 1; + nsplit2m1 = nsplit2  1; + xsplit2[nsplit2m1  k] = xintersect[0]; + ysplit2[nsplit2m1  k] = yintersect[0]; + ifsplit2[nsplit2m1  k] = 1; + kstart = k + 1; + kk = kkstart1; /* No wrap checks on kk index below because * it must always be in valid range (since * polygon 1 traversed only once). */ for ( k = kstart; k < range1 + 1; k++ ) {  xsplit1[k] = x1[kk];  ysplit1[k] = y1[kk];  ifsplit1[k] = 2;  xsplit2[k] = x1[kk];  ysplit2[k] = y1[kk++];  ifsplit2[k] = 2; + xsplit1[k] = x1[kk]; + ysplit1[k] = y1[kk]; + ifsplit1[k] = 2; + xsplit2[nsplit2m1  k] = x1[kk]; + ysplit2[nsplit2m1  k] = y1[kk++]; + ifsplit2[nsplit2m1  k] = 2; }  xsplit1[k] = xintersect[1];  ysplit1[k] = yintersect[1];  ifsplit1[k] = 1;  xsplit2[k] = xintersect[1];  ysplit2[k] = yintersect[1];  ifsplit2[k] = 1; + xsplit1[k] = xintersect[1]; + ysplit1[k] = yintersect[1]; + ifsplit1[k] = 1; + xsplit2[nsplit2m1  k] = xintersect[1]; + ysplit2[nsplit2m1  k] = yintersect[1]; + ifsplit2[nsplit2m1  k] = 1; /* Finish off collecting split1 using ascending kk * values. */ @@ 1608,6 +1667,11 @@ kk = n2; } + /* N.B. the positive orientation of split1 is + * preserved since the index order is the same + * as that of polygon 2, and by assumption + * that polygon and polygon 1 have identical + * positive orientations. */ fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, x1, y1, i1start, n1, @@ 1621,13 +1685,18 @@ kk = kkstart22; for ( k = kstart; k < nsplit2; k++ ) {  xsplit2[k] = x2[kk];  ysplit2[k] = y2[kk];  ifsplit2[k] = if2[kk]; + xsplit2[nsplit2m1  k] = x2[kk]; + ysplit2[nsplit2m1  k] = y2[kk]; + ifsplit2[nsplit2m1  k] = if2[kk]; if ( kk < 0 ) kk += n2; } + /* N.B. the positive orientation of split2 is + * preserved since the index order is the same + * as that of polygon 2, and by assumption + * that polygon and polygon 1 have identical + * positive orientations. */ fill_intersection_polygon( recursion_depth + 1, ifextrapolygon, fill, x1, y1, i1start, n1, @@ 1793,9 +1862,9 @@ { /* Two line segments are parallel */ status = status  PL_PARALLEL;  /* Choice of intersect is arbitrary in this case. Choose A1 if  * that lies near or in B. Otherwise, choose A2 if that lies near  * or in B. Otherwise, choose the average point. */ + /* Choice of intersect is arbitrary in this case. Choose A1, A2, + * B1, or B2 (in that order) if any of them lie inside or near + * the other line segment. Otherwise, choose the average point. */ if (( BETW_NBCC( xA1, xB1, xB2 ) && BETW_NBCC( yA1, yB1, yB2 ))) { fxintersect = xA1; @@ 1806,6 +1875,16 @@ fxintersect = xA2; fyintersect = yA2; } + else if (( BETW_NBCC( xB1, xA1, xA2 ) && BETW_NBCC( yB1, yA1, yA2 ))) + { + fxintersect = xB1; + fyintersect = yB1; + } + else if (( BETW_NBCC( xB2, xA1, xA2 ) && BETW_NBCC( yB2, yA1, yA2 ))) + { + fxintersect = xB2; + fyintersect = yB2; + } else { fxintersect = 0.25 * ( xA1 + xA2 + xB1 + xB2 ); @@ 1849,3 +1928,32 @@ return status; } + +/* Decide if polygon has a positive orientation or not. + * See http://en.wikipedia.org/wiki/Curve_orientation for details + * of this simple determinate method. */ +int +positive_orientation( PLINT n, const PLINT *x, const PLINT *y ) +{ + PLFLT xa, ya, xb, yb, xc, yc, det; + if ( n < 3 ) + { + plwarn( "positive_orientation: internal logic error, n < 3" ); + return 0; + } + /* Use floating point to avoid integer overflows. */ + xa = x[0]; + xb = x[1]; + xc = x[2]; + ya = y[0]; + yb = y[1]; + yc = y[2]; + det = ( xa * yb + xb * yc + xc * ya )  ( xa * yc + xb * ya + xc * yb ); + if ( det == 0. ) + { + plwarn( "positive_orientation: internal logic error, det == 0." ); + return 0; + } + else + return det > 0.; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091222 19:15:29

Revision: 10728 http://plplot.svn.sourceforge.net/plplot/?rev=10728&view=rev Author: airwin Date: 20091222 19:15:21 +0000 (Tue, 22 Dec 2009) Log Message:  Upgrade positive_orientation logic from a simple rule which only works consistently for convex polygons to a general rule which works for all nonintersecting polygons whether convex or not. By chance, this DUSE_FILL_INTERSECTION_POLYGON=ON bug fix makes no difference in the results for example 25, but it will make a difference in the general case. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091222 02:31:06 UTC (rev 10727) +++ trunk/src/plfill.c 20091222 19:15:21 UTC (rev 10728) @@ 1930,30 +1930,40 @@ } /* Decide if polygon has a positive orientation or not.  * See http://en.wikipedia.org/wiki/Curve_orientation for details  * of this simple determinate method. */ + * Note the simple algorithm given in + * http://en.wikipedia.org/wiki/Curve_orientation is incorrect for + * nonconvex polygons. Instead, for the general nonintersecting case + * use the polygon area method given by + * http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/ where + * you project each edge of the polygon down to the X axis and take the + * area of the enclosed trapezoid. The trapezoid areas outside the + * polygon are completely cancelled if you sum over all edges. Furthermore, + * the sum of the trapezoid areas have terms which are zero when calculated + * with the telescoping rule so the final formula is quite simple. */ int positive_orientation( PLINT n, const PLINT *x, const PLINT *y ) {  PLFLT xa, ya, xb, yb, xc, yc, det; + PLINT i, im1; + /* Use PLFLT for all calculations to avoid integer overflows. Also, + * the normal PLFLT has 64bits which means you get exact integer + * arithmetic well beyond the normal integer overflow limits. */ + PLFLT twice_area = 0.; if ( n < 3 ) { plwarn( "positive_orientation: internal logic error, n < 3" ); return 0; }  /* Use floating point to avoid integer overflows. */  xa = x[0];  xb = x[1];  xc = x[2];  ya = y[0];  yb = y[1];  yc = y[2];  det = ( xa * yb + xb * yc + xc * ya )  ( xa * yc + xb * ya + xc * yb );  if ( det == 0. ) + im1 = n  1; + for ( i = 0; i < n; i++ ) {  plwarn( "positive_orientation: internal logic error, det == 0." ); + twice_area += (PLFLT) x[im1] * (PLFLT) y[i]  (PLFLT) x[i] * (PLFLT) y[im1]; + im1 = i; + } + if ( twice_area == 0. ) + { + plwarn( "positive_orientation: internal logic error, twice_area == 0." ); return 0; } else  return det > 0.; + return twice_area > 0.; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <arjenmarkus@us...>  20091224 09:19:37

Revision: 10729 http://plplot.svn.sourceforge.net/plplot/?rev=10729&view=rev Author: arjenmarkus Date: 20091224 09:19:26 +0000 (Thu, 24 Dec 2009) Log Message:  Needed to copy the check on npts and draw() because the (ancient) MS Visual C 6.0 compiler complained. Copied the check to after the declarations. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091222 19:15:21 UTC (rev 10728) +++ trunk/src/plfill.c 20091224 09:19:26 UTC (rev 10729) @@ 456,8 +456,6 @@ PLINT xmin, PLINT xmax, PLINT ymin, PLINT ymax, void ( *draw )( short *, short *, PLINT )) {  /* Must have at least 3 points and draw() specified */  if ( npts < 3  !draw ) return; #ifdef USE_FILL_INTERSECTION_POLYGON PLINT *x10, *y10, *x1, *y1, i1start = 0, i, im1, n1, n1m1; @@ 465,6 +463,10 @@ PLINT y2[4] = { ymin, ymin, ymax, ymax }; PLINT if2[4] = { 0, 0, 0, 0 }; PLINT n2 = 4; + + /* Must have at least 3 points and draw() specified */ + if ( npts < 3  !draw ) return; + if (( x10 = (PLINT *) malloc( npts * sizeof ( PLINT ))) == NULL ) { plexit( "plP_plfclp: Insufficient memory" ); @@ 549,6 +551,8 @@ int inside_rb; int inside_ru; + /* Must have at least 3 points and draw() specified */ + if ( npts < 3  !draw ) return; if ( npts < PL_MAXPOLY ) { This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091227 19:27:19

Revision: 10730 http://plplot.svn.sourceforge.net/plplot/?rev=10730&view=rev Author: airwin Date: 20091227 19:27:12 +0000 (Sun, 27 Dec 2009) Log Message:  Implement PL_NEAR_PARALLEL crossing status. Implement fill_status to keep track of whether a given split of polygon 2 is known to have filled areas or known to not have filled areas. At the final stage where there are not crossings left between the polygons, this status helps to quickly decide whether to fill or not. Prepare calls to fill_intersection_polygon so that polygon 1 always starts outside polygon 2. That and the positive orientation requirement considerably simplify the logic of fill_intersection_polygon. Implement number_crossings to make sure that the crossings are determined in the correct order for a give segment of polygon 1. N.B. These changes build, but still produce incorrect fill results for some of the subpages of example 25. So this commit is in the nature of a code drop to keep track of the many changes that have been done with the DUSE_FILL_INTERSECTION_POLYGON=ON approach rather than a final result. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091224 09:19:26 UTC (rev 10729) +++ trunk/src/plfill.c 20091227 19:27:12 UTC (rev 10730) @@ 36,8 +36,8 @@ /* Status codes ORed together in the return from notcrossed. * PL_NOT_CROSSED is set whenever the two line segments definitely  * (i.e., not near the ends or completely apart) do not cross each  * other. + * (i.e., intersection not near the ends or completely apart) + * do not cross each other. * * PL_NEAR_A1 is set if the intersect is near (+/ PL_NBCC) the first * end of line segment A. @@ 51,16 +51,22 @@ * PL_NEAR_B2 is set if the intersect is near (+/ PL_NBCC) the second * end of line segment B. *  * PL_PARALLEL is set if the two line segments are parallel with each other. + * PL_NEAR_PARALLEL is set if the two line segments are nearly parallel + * with each other, i.e., a change in one of the coordinates of up to + * (+/ PL_NBCC) would render them exactly parallel. + * + * PL_PARALLEL is set if the two line segments are exactly parallel + * with each other. */ enum PL_CrossedStatus {  PL_NOT_CROSSED = 0x1,  PL_NEAR_A1 = 0x2,  PL_NEAR_A2 = 0x4,  PL_NEAR_B1 = 0x8,  PL_NEAR_B2 = 0x10,  PL_PARALLEL = 0x20 + PL_NOT_CROSSED = 0x1, + PL_NEAR_A1 = 0x2, + PL_NEAR_A2 = 0x4, + PL_NEAR_B1 = 0x8, + PL_NEAR_B2 = 0x10, + PL_NEAR_PARALLEL = 0x20, + PL_PARALLEL = 0x40 }; struct point @@ 91,6 +97,7 @@ static void fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, + PLINT fill_status, void ( *fill )( short *, short *, PLINT ), const PLINT *x1, const PLINT *y1, PLINT i1start, PLINT n1, @@ 105,6 +112,12 @@ static int positive_orientation( PLINT n, const PLINT *x, const PLINT *y ); +static int +number_crossings( PLINT *xcross, PLINT *ycross, PLINT *i2cross, PLINT ncross, + PLINT i1, PLINT n1, const PLINT *x1, const PLINT *y1, + PLINT n2, const PLINT *x2, const PLINT *y2 ); + + /**\ * void plfill() * @@ 456,9 +469,9 @@ PLINT xmin, PLINT xmax, PLINT ymin, PLINT ymax, void ( *draw )( short *, short *, PLINT )) {  #ifdef USE_FILL_INTERSECTION_POLYGON  PLINT *x10, *y10, *x1, *y1, i1start = 0, i, im1, n1, n1m1; + PLINT *x10, *y10, *x1, *y1, *if1, i1start = 0, i, im1, n1, n1m1, + ifnotpointinpolygon; PLINT x2[4] = { xmin, xmax, xmax, xmin }; PLINT y2[4] = { ymin, ymin, ymax, ymax }; PLINT if2[4] = { 0, 0, 0, 0 }; @@ 528,7 +541,30 @@ free( x10 ); free( y10 ); }  fill_intersection_polygon( 0, 0, draw, x1, y1, i1start, n1, x2, y2, if2, n2 ); + + /* Insure that the first vertex of polygon 1 (starting with n1  + * 1) that is not on the border of polygon 2 is definitely outside + * polygon 2. */ + im1 = n1  1; + for ( i = 0; i < n1; i++ ) + { + if (( ifnotpointinpolygon = + notpointinpolygon( n2, x2, y2, x1[im1], y1[im1] )) != 1 ) + break; + im1 = i; + } + + if ( ifnotpointinpolygon ) + fill_intersection_polygon( 0, 0, 0, draw, x1, y1, i1start, n1, x2, y2, if2, n2 ); + else + { + if (( if1 = (PLINT *) calloc( n1, sizeof ( PLINT ))) == NULL ) + { + plexit( "plP_plfclp: Insufficient memory" ); + } + fill_intersection_polygon( 0, 0, 0, draw, x2, y2, i1start, n2, x1, y1, if1, n1 ); + free( if1 ); + } free( x1 ); free( y1 ); return; @@ 1297,15 +1333,17 @@ #define MAX_RECURSION_DEPTH 10 /* Fill intersection of two simple (not selfintersecting) polygons  * that both have a positive orientation (see +/* Fill intersection of two simple polygons that do no selfintersect, + * and which have no duplicate vertices or two consecutive edges that + * are parallel. A further requirement is that both polygons have a + * positive orientation (see * http://en.wikipedia.org/wiki/Curve_orientation). That is, as you * traverse the boundary in index order, the inside area of the  * polygon is always on the left. This requirement simplifies the  * logic of fill_instersection_polygon. N.B. it is the calling  * routine's responsibility to insure the two polygons do not have  * duplicate points, are not selfintersecting, and have positive  * orientation. + * polygon is always on the left. Finally, the first vertex of + * polygon 1 (starting with n1 1) that is not near the border of + * polygon 2 must be outside polygon 2. N.B. it is the calling + * routine's responsibility to insure all those requirements are + * satisfied. * * Two polygons that do not self intersect must have an even number of * edge crossings between them. (ignoring vertex intersections which @@ 1322,20 +1360,41 @@ * ifextrapolygon false, and the second polygon set to an additional * polygon defined by the stream (not yet implemented). */ +/* arguments to intersection_polygon: + * recursion_depth is just what it says. + * ifextrapolygon used to decide whether to use extra polygon from the stream. + *fill is the fill routine. + *x1, *y1, n1 define the polygon 1 vertices. + * i1start is the first polygon 1 index to look at (because all the previous + * ones have been completely processed). + *x2, *y2, *if2, n2 define the polygon 2 vertices plus a status indicator + * for each vertex which is 1 for a previous crossing and 2 for a polygon + * 1 vertex. + * fill_status is 1 when polygons 1 and 2 _must_ include some joint + * filled area and is 1 when polygons 1 and 2 _must_ include some + * unfilled area. fill_status of +/ 1 is determined from the + * orientations of polygon 1 and 2 from the next higher recursion + * level and how those two are combined to form the polygon 2 + * split at this recursion level. fill_status = 0 occurs (at + * recursion level 0) for polygons 1 and 2 that are independent of + * each other. */ + void fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon, + PLINT fill_status, void ( *fill )( short *, short *, PLINT ), const PLINT *x1, const PLINT *y1, PLINT i1start, PLINT n1, const PLINT *x2, const PLINT *y2, const PLINT *if2, PLINT n2 ) {  PLINT i1, i1m1, i1wrap, i1wraplast,  i2, i2m1, i2wrap, i2wraplast, + PLINT i1, i1m1, i1start_new, + i2, i2m1, kk, kkstart1, kkstart21, kkstart22,  k, kstart, range1, range21, range22, ncrossed, + k, kstart, range1, range21, range22, ncrossed, ncrossed_change, nsplit1, nsplit2, nsplit2m1;  PLINT xintersect[2], yintersect[2], ifnotcrossed; + PLINT xintersect[2], yintersect[2], i1intersect[2], + i2intersect[2], ifnotcrossed; PLINT *xsplit1, *ysplit1, *ifsplit1, *xsplit2, *ysplit2, *ifsplit2; PLINT ifill, nfill = 0, @@ 1368,6 +1427,7 @@ return; } + /* Check that there are no duplicate vertices. */ i1m1 = i1start  1; if ( i1m1 < 0 ) i1m1 = n1  1; @@ 1455,268 +1515,265 @@ * split, and in that case we move on to the end phase below. */ ncrossed = 0;  i1wrap = i1start  1;  if ( i1wrap < 0 )  i1m1 = n1  1;  else  i1m1 = i1wrap; + i1m1 = i1start  1; + if ( i1m1 < 0 ) + i1m1 += n1; for ( i1 = i1start; i1 < n1; i1++ ) {  i2m1 = n2  1;  i2wrap = 1;  for ( i2 = 0; i2 < n2; i2++ ) + i2m1 = n2  1; + ncrossed_change = number_crossings( + &xintersect[ncrossed], &yintersect[ncrossed], + &i2intersect[ncrossed], 2  ncrossed, + i1, n1, x1, y1, n2, x2, y2 ); + if ( ncrossed_change > 0 ) {  if ( !( if2[i2] && if2[i2m1] )) + i1intersect[ncrossed] = i1; + if ( ncrossed_change == 2 ) ; + i1intersect[1] = i1; + + ncrossed += ncrossed_change; + if ( ncrossed == 2 ) {  ifnotcrossed = notcrossed(  &xintersect[ncrossed], &yintersect[ncrossed],  x1[i1m1], y1[i1m1], x1[i1], y1[i1],  x2[i2m1], y2[i2m1], x2[i2], y2[i2] );  /* Use only definite crossing case. */  if ( !ifnotcrossed ) + /* Have discovered the first two crossings for + * polygon 1 at i1 = i1start or above. */ + + /* Calculate polygon 2 index range in split 1 (the + * split that proceeds beyond the second intersect with + * ascending i2 values). */ + range21 = i2intersect[0]  i2intersect[1]; + if ( range21 < 0 ) + range21 += n2; + /* i2 intersect values range between 0 and n2  1 so + * the the smallest untransformed range21 value is n2 + + 1, and the largest untransformed range21 value is + * n2  1. This means the smallest transformed range21 + * value is +1. Also, for the untransformed range21 + * value of 0, go the "long way" around with + * transformed range21 = n2 so that split 1 is + * guaranteed to have the correct fill_status == +1. */ + + /* N.B. always avoid special code for range21 == 0 + * case because it is no longer necessary since there + * is a guarantee that the first vertex of polygon 1 + * (starting with n1  1) that is not near the border + * of polygon 2 will be outside polygon 2. */ + if ( 0 && range21 == 0 ) {  if ( ncrossed == 0 ) + /* For this special case the above ascii art + * does not apply, and we must use the + * following diagram instead: + * + *  1 1 ... + * + * + * + * 2???2 X X 2???2 + * + * 1 1 + * + * + * N.B. no valid split of polygon 2 is + * possible with this diagram. However, + * this diagrem can be classified the same + * as the first diagram, but with the roles + * of polygon 1 and polygon 2 swapped so + * that polygon 1 is the one that gets split + * instead of polygon 2. In fact, swapping + * those two polygons (and starting the new + * polygon 1 just at i2intersect[1] of the old + * polygon 2 to insure the present two + * intersections are immediately (at least for the + * first new polygon segment 1) found and dealt + * with to eliminate the chance of an infinite + * "swap" loop) is the only clean way to deal with + * this situation. */ + if (( xsplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL ) {  i1wraplast = i1wrap;  i2wraplast = i2wrap;   ncrossed++; + plexit( "fill_intersection_polygon: Insufficient memory" ); }  else + if (( ysplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL ) {  /* Have discovered the first two intersections for  * polygon 1 at i1 = i1start or above. */ + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ifsplit1 = (PLINT *) calloc( n1, sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + kk = i2intersect[1]; + for ( i2 = 0; i2 < n2; i2++ ) + { + xsplit2[kk] = x2[i2]; + ysplit2[kk++] = y2[i2]; + if ( kk == n2 ) + kk = 0; + } + /* N.B. the positive orientation of split2 + * is preserved since the index order is + * the same as that of polygon 2, and by + * assumption that polygon and polygon 1 + * have identical positive + * orientations. fill_status does not matter + * here because a split is guaranteed at the + * next recursion level. For definiteness set + * it to zero. */ + fill_intersection_polygon( + recursion_depth + 1, ifextrapolygon, + 0, fill, + xsplit2, ysplit2, 0, n2, + x1, y1, ifsplit1, n1 ); + free( xsplit2 ); + free( ysplit2 ); + free( ifsplit1 ); + return; + } + /* New i1start is the same as the current i1 (just + * in case there are more crossings to find between + * i1m1 and i1.) */ + i1start_new = i1;  /* Calculate polygon 2 index range in split 1  * (the split that normally proceeds beyond  * the second intersect with ascending i2  * values). */  range21 = i2wraplast  i2wrap;  if ( range21 < 0 )  range21 += n2;  /* Wrap values range between 1 and n2  2 so  * the the smallest untransformed range21  * value is n2 + 1, the smallest transformed  * range21 value is 0 (only obtained if  * i2wraplast == i2wrap), and the largest  * transformed range21 value is n2  1. */  if ( range21 == 0 )  {  /* For this special case the above ascii art  * does not apply, and we must use the  * following diagram instead:  *  *  1 1 ...  *  *  *  * 2???2 X X 2???2  *  * 1 1  *  *  * N.B. no valid split of polygon 2 is  * possible with this diagram. However,  * this diagrem can be classified the same  * as the first diagram, but with the roles  * of polygon 1 and polygon 2 swapped so  * that polygon 1 is the one that gets split  * instead of polygon 2. In fact, swapping  * those two polygons (and starting the new  * polygon 1 just at i2 of the old polygon 2  * to insure the present two intersections  * are immediately (at least for the first  * new polygon segment 1) found and dealt  * with to eliminate the chance of an  * infinite "swap" loop) is the only clean  * way to deal with this situation. */  if (( xsplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ysplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ifsplit1 = (PLINT *) calloc( n1, sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  kk = i2;  for ( i2 = 0; i2 < n2; i2++ )  {  xsplit2[kk] = x2[i2];  ysplit2[kk++] = y2[i2];  if ( kk == n2 )  kk = 0;  }  /* N.B. the positive orientation of split2  * is preserved since the index order is  * the same as that of polygon 2, and by  * assumption that polygon and polygon 1  * have identical positive  * orientations. */  fill_intersection_polygon(  recursion_depth + 1, ifextrapolygon, fill,  xsplit2, ysplit2, 0, n2,  x1, y1, ifsplit1, n1 );  free( xsplit2 );  free( ysplit2 );  free( ifsplit1 );  return;  }  /* New i1start is the largest nonnegative polygon 1  * index below the last detected intersect. */  i1start = MAX( i1wrap, 0 ); + /* Split polygon 2 at the boundary consisting of + * first intersection, intervening (if any) range1 + * polygon 1 points and second intersection. */ + /* range1 must always be nonnegative because i1 + * range only traversed once. */ + range1 = i1intersect[1]  i1intersect[0]; + /* Polygon 2 intersects could be anywhere (since + * i2 range repeated until get an intersect). + * Divide polygon 2 into two polygons with a + * common boundary consisting of the first intersect, + * range1 points from polygon 1 starting at index + * kkstart1 of polygon 1, and the second intersect. */ + kkstart1 = i1intersect[0];  /* Split polygon 2 at the boundary consisting of  * first intersection, intervening (if any) range1  * polygon 1 points and second intersection. */  /* range1 must always be nonnegative because i1  * range only traversed once. */  range1 = i1wrap  i1wraplast;  /* Polygon 2 intersects could be anywhere (since  * i2 range repeated until get an intersect).  * Divide polygon 2 into two polygons with a  * common boundary consisting of the first intersect,  * range1 points from polygon1 starting at index  * kkstart1 of polygon 1, and the second intersect. */  kkstart1 = i1wraplast + 1; + kkstart21 = i2intersect[1];  kkstart21 = i2; + /* Split 2 of polygon 2 consists of the + * boundary + range22 (= n2  range21) points + * between kkstart22 (= i2intersect[1]1) and i2intersect[0] in + * descending order of polygon 2 indices. */ + range22 = n2  range21; + kkstart22 = i2intersect[1]  1; + if ( kkstart22 < 0 ) + kkstart22 += n2; + nsplit1 = 2 + range1 + range21; + nsplit2 = 2 + range1 + range22;  /* Split 2 of polygon 2 consists of the  * boundary + range22 (= n2  range21) points  * between kkstart22 (= i2wrap) and i2wraplast + 1 in  * descending order of polygon 2 indices. */  range22 = n2  range21;  kkstart22 = i2wrap;  if ( kkstart22 == 1 )  kkstart22 = i2m1;  nsplit1 = 2 + range1 + range21;  nsplit2 = 2 + range1 + range22; + if (( xsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ysplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ifsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + }  if (( xsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ysplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ifsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  } + if (( xsplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ysplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + if (( ifsplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ))) == NULL ) + { + plexit( "fill_intersection_polygon: Insufficient memory" ); + } + /* Common boundary between split1 and split2. */ + /* N.B. Although basic index arithmetic for + * split 2 is done in negative orientation + * order because the index is decrementing + * relative to the index of split 2, actually + * store results in reverse order to preserve + * the positive orientation that by assumption + * both polygon 1 and 2 have. */ + k = 0; + xsplit1[k] = xintersect[0]; + ysplit1[k] = yintersect[0]; + ifsplit1[k] = 1; + nsplit2m1 = nsplit2  1; + xsplit2[nsplit2m1  k] = xintersect[0]; + ysplit2[nsplit2m1  k] = yintersect[0]; + ifsplit2[nsplit2m1  k] = 1; + kstart = k + 1; + kk = kkstart1; + /* No wrap checks on kk index below because + * it must always be in valid range (since + * polygon 1 traversed only once). */ + for ( k = kstart; k < range1 + 1; k++ ) + { + xsplit1[k] = x1[kk]; + ysplit1[k] = y1[kk]; + ifsplit1[k] = 2; + xsplit2[nsplit2m1  k] = x1[kk]; + ysplit2[nsplit2m1  k] = y1[kk++]; + ifsplit2[nsplit2m1  k] = 2; + } + xsplit1[k] = xintersect[1]; + ysplit1[k] = yintersect[1]; + ifsplit1[k] = 1; + xsplit2[nsplit2m1  k] = xintersect[1]; + ysplit2[nsplit2m1  k] = yintersect[1]; + ifsplit2[nsplit2m1  k] = 1;  if (( xsplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ysplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ifsplit2 = (PLINT *) malloc( nsplit2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  /* Common boundary between split1 and split2. */  /* N.B. Although basic index arithmetic for  * split 2 is done in negative orientation  * order because the index is decrementing  * relative to the index of split 2, actually  * store results in reverse order to preserve  * the positive orientation that by assumption  * both polygon 1 and 2 have. */  k = 0;  xsplit1[k] = xintersect[0];  ysplit1[k] = yintersect[0];  ifsplit1[k] = 1;  nsplit2m1 = nsplit2  1;  xsplit2[nsplit2m1  k] = xintersect[0];  ysplit2[nsplit2m1  k] = yintersect[0];  ifsplit2[nsplit2m1  k] = 1;  kstart = k + 1;  kk = kkstart1;  /* No wrap checks on kk index below because  * it must always be in valid range (since  * polygon 1 traversed only once). */  for ( k = kstart; k < range1 + 1; k++ )  {  xsplit1[k] = x1[kk];  ysplit1[k] = y1[kk];  ifsplit1[k] = 2;  xsplit2[nsplit2m1  k] = x1[kk];  ysplit2[nsplit2m1  k] = y1[kk++];  ifsplit2[nsplit2m1  k] = 2;  }  xsplit1[k] = xintersect[1];  ysplit1[k] = yintersect[1];  ifsplit1[k] = 1;  xsplit2[nsplit2m1  k] = xintersect[1];  ysplit2[nsplit2m1  k] = yintersect[1];  ifsplit2[nsplit2m1  k] = 1; + /* Finish off collecting split1 using ascending kk + * values. */ + kstart = k + 1; + kk = kkstart21; + for ( k = kstart; k < nsplit1; k++ ) + { + xsplit1[k] = x2[kk]; + ysplit1[k] = y2[kk]; + ifsplit1[k] = if2[kk++]; + if ( kk >= n2 ) + kk = n2; + }  /* Finish off collecting split1 using ascending kk  * values. */  kstart = k + 1;  kk = kkstart21;  for ( k = kstart; k < nsplit1; k++ )  {  xsplit1[k] = x2[kk];  ysplit1[k] = y2[kk];  ifsplit1[k] = if2[kk++];  if ( kk >= n2 )  kk = n2;  } + /* N.B. the positive orientation of split1 is + * preserved since the index order is the same + * as that of polygon 2, and by assumption + * that polygon and polygon 1 have identical + * positive orientations. */ + fill_intersection_polygon( + recursion_depth + 1, ifextrapolygon, 1, fill, + x1, y1, i1start_new, n1, + xsplit1, ysplit1, ifsplit1, nsplit1 ); + free( xsplit1 ); + free( ysplit1 ); + free( ifsplit1 );  /* N.B. the positive orientation of split1 is  * preserved since the index order is the same  * as that of polygon 2, and by assumption  * that polygon and polygon 1 have identical  * positive orientations. */  fill_intersection_polygon(  recursion_depth + 1, ifextrapolygon, fill,  x1, y1, i1start, n1,  xsplit1, ysplit1, ifsplit1, nsplit1 );  free( xsplit1 );  free( ysplit1 );  free( ifsplit1 ); + /* Finish off collecting split2 using descending kk + * values. */ + kk = kkstart22; + for ( k = kstart; k < nsplit2; k++ ) + { + xsplit2[nsplit2m1  k] = x2[kk]; + ysplit2[nsplit2m1  k] = y2[kk]; + ifsplit2[nsplit2m1  k] = if2[kk]; + if ( kk < 0 ) + kk += n2; + }  /* Finish off collecting split2 using descending kk  * values. */  kk = kkstart22;  for ( k = kstart; k < nsplit2; k++ )  {  xsplit2[nsplit2m1  k] = x2[kk];  ysplit2[nsplit2m1  k] = y2[kk];  ifsplit2[nsplit2m1  k] = if2[kk];  if ( kk < 0 )  kk += n2;  }   /* N.B. the positive orientation of split2 is  * preserved since the index order is the same  * as that of polygon 2, and by assumption  * that polygon and polygon 1 have identical  * positive orientations. */  fill_intersection_polygon(  recursion_depth + 1, ifextrapolygon, fill,  x1, y1, i1start, n1,  xsplit2, ysplit2, ifsplit2, nsplit2 );  free( xsplit2 );  free( ysplit2 );  free( ifsplit2 );  return;  }  } + /* N.B. the positive orientation of split2 is + * preserved since the index order is the same + * as that of polygon 2, and by assumption + * that polygon and polygon 1 have identical + * positive orientations. */ + fill_intersection_polygon( + recursion_depth + 1, ifextrapolygon, 1, fill, + x1, y1, i1start_new, n1, + xsplit2, ysplit2, ifsplit2, nsplit2 ); + free( xsplit2 ); + free( ysplit2 ); + free( ifsplit2 ); + return; }  i2m1 = i2;  i2wrap = i2m1; }  i1m1 = i1;  i1wrap = i1m1; + i1m1 = i1; } if ( ncrossed != 0 ) @@ 1725,73 +1782,91 @@ return; }  /* This end phase is reached only if no intersections are found.  * For this case, there are 5 valid and 1 invalid combination of  * the following conditions:the intersection of polygon 2 and 1, must be  * either of them (in which case fill with the inner one), or neither  * of them (in which case don't fill at all). */ + /* This end phase is reached only if no crossings are found. */  /* Classify polygon 1 by looking for first vertex in polygon 1  * that is definitely inside or outside polygon 2. */  for ( i1 = 0; i1 < n1; i1++ ) + /* If a fill_status of +/ 1 is known, use that to fill or not since + +1 corresponds to all of polygon 2 inside polygon 1 and 1 + * corresponds to none of polygon 2 inside polygon 1. */ + if ( fill_status == 1 ) + return; + else if ( fill_status == 1 ) {  if (( ifnotpolygon1inpolygon2 =  notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] )) != 1 )  break; + nfill = n2; + xfiller = x2; + yfiller = y2; }   /* Classify polygon 2 by looking for first vertex in polygon 2  * that is definitely inside or outside polygon 1. */  ifnotpolygon2inpolygon1 = 1;  for ( i2 = 0; i2 < n2; i2++ ) + else if ( fill_status == 0 ) + //else if ( 1 ) {  /* Do not bother checking vertices already known to be on the  * boundary with polygon 1. */  if ( !if2[i2] && ( ifnotpolygon2inpolygon1 =  notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] )) != 1 )  break;  } + if ( recursion_depth != 0 ) + { + plwarn( "fill_intersection_polygon: Internal error; fill_status == 0 for recursion_depth > 0" ); + return; + } + /* For this case (recursion level 0) the two polygons are + * completely independent with no crossings between them or + * edges constructed from one another. + * + * The intersection of polygon 2 and 1, must be either of them (in + * which case fill with the inner one), or neither of them (in + * which case don't fill at all). */  if ( ifnotpolygon2inpolygon1 == 0 )  {  /* Polygon 2 definitely inside polygon 1. */  if ( ifnotpolygon1inpolygon2 != 0 ) + /* Classify polygon 1 by looking for first vertex in polygon 1 + * that is definitely inside or outside polygon 2. */ + for ( i1 = 0; i1 < n1; i1++ ) {  nfill = n2;  xfiller = x2;  yfiller = y2; + if (( ifnotpolygon1inpolygon2 = + notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] )) != 1 ) + break; }  else + + /* Classify polygon 2 by looking for first vertex in polygon 2 + * that is definitely inside or outside polygon 1. */ + ifnotpolygon2inpolygon1 = 1; + for ( i2 = 0; i2 < n2; i2++ ) + { + /* Do not bother checking vertices already known to be on the + * boundary with polygon 1. */ + if ( !if2[i2] && ( ifnotpolygon2inpolygon1 = + notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] )) != 1 ) + break; + } + + if ( ifnotpolygon2inpolygon1 == 0 && ifnotpolygon1inpolygon2 == 0 ) plwarn( "fill_intersection_polygon: Internal error; no intersections found but each polygon definitely inside the other!" );  }  else if ( ifnotpolygon2inpolygon1 == 1 )  {  /* Polygon 2 vertices near polygon 1 border. */  if ( ifnotpolygon1inpolygon2 != 0 ) + else if ( ifnotpolygon2inpolygon1 == 2 && ifnotpolygon1inpolygon2 == 2 ) + /* The polygons do not intersect each other so do not fill in this + * case. */ + return; + else if ( ifnotpolygon2inpolygon1 == 0 ) {  /* Polygon 1 vertices near polygon 2 border or definitely outside it. */ + /* Polygon 2 definitely inside polygon 1. */ nfill = n2; xfiller = x2; yfiller = y2; }  else + else if ( ifnotpolygon1inpolygon2 == 0 ) {  /* Polygon 1 vertices definitely inside polygon 2 border. */ + /* Polygon 1 definitely inside polygon 2. */ nfill = n1; xfiller = x1; yfiller = y1; }  }  else  {  /* Polygon 2 vertices definitely outside polygon 1 border. */  if ( ifnotpolygon1inpolygon2 != 2 ) + else if ( ifnotpolygon2inpolygon1 == 1 && ifnotpolygon1inpolygon2 == 1 ) {  /* Polygon 1 vertices definitely inside polygon 2 border. */  nfill = n1;  xfiller = x1;  yfiller = y1; + /* Polygon 2 vertices near polygon 1 border and vice versa which + * implies the polygons are identical. */ + nfill = n2; + xfiller = x2; + yfiller = y2; } + else + { + /* Polygon 1 inscribed in polygon 2 or vice versa. This is normally + * unlikely for two independent polygons so the implementation is + * ToDo. */ + plwarn( "fill_intersection_polygon: inscribed polygons are still ToDo" ); + } } if ( nfill > 0 ) @@ 1832,7 +1907,7 @@ PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2, PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 ) {  PLFLT factor, fxintersect, fyintersect; + PLFLT factor, factor_NBCC, fxintersect, fyintersect; /* These variables are PLFLT not for precision, but to * avoid integer overflows if they were typed as PLINT. */ PLFLT xA2A1, yA2A1, xB2B1, yB2B1; @@ 1861,11 +1936,14 @@ xB2B1 = xB2  xB1; yB2B1 = yB2  yB1;  factor = xA2A1 * yB2B1  yA2A1 * xB2B1;  if ( fabs( factor ) == 0. ) + factor = xA2A1 * yB2B1  yA2A1 * xB2B1; + factor_NBCC = PL_NBCC * ( fabs( xA2A1 ) + fabs( yB2B1 ) + fabs( yA2A1 ) + fabs( xB2B1 )); + if ( fabs( factor ) <= factor_NBCC ) {  /* Two line segments are parallel */  status = status  PL_PARALLEL; + if ( fabs( factor ) > 0. ) + status = status  PL_NEAR_PARALLEL; + else + status = status  PL_PARALLEL; /* Choice of intersect is arbitrary in this case. Choose A1, A2, * B1, or B2 (in that order) if any of them lie inside or near * the other line segment. Otherwise, choose the average point. */ @@ 1971,3 +2049,103 @@ else return twice_area > 0.; } + +/* For the line with endpoints which are the (i11)th, and i1th + * vertices of polygon 1 with polygon 2 find all definite crossings + * with polygon 1. (The full polygon 1 information is needed only to + * help sort out (NOT IMPLEMENTED YET) ambiguous crossings near a + * vertex of polygon 1.) Sort those definite crossings in ascending + * order along the line between the (i11)th and i1th vertices of + * polygon 1, and return the first ncross (= 1 or = 2) crossings in the + * xcross, ycross, and i2cross arrays. Return a zero or positive + * status code of the actual number of crossings retained up to the + * maximum allowed value of ncross. If some error occurred, return a + * negative status code. */ + +int +number_crossings( PLINT *xcross, PLINT *ycross, PLINT *i2cross, PLINT ncross, + PLINT i1, PLINT n1, const PLINT *x1, const PLINT *y1, + PLINT n2, const PLINT *x2, const PLINT *y2 ) +{ + int i1m1, i2, i2m1, ifnotcrossed; + int ifxsort, ifascend, count_crossings = 0, status = 0; + PLINT xintersect, yintersect; + + i1m1 = i1  1; + if ( i1m1 < 0 ) + i1m1 += n1; + if ( !( ncross == 1  ncross == 2 )  + ( x1[i1m1] == x1[i1] && y1[i1m1] == y1[i1] )  n1 < 2  n2 < 2 ) + { + plwarn( "findcrossings: invalid call" ); + return 1; + } + + ifxsort = abs( x1[i1]  x1[i1m1] ) > abs( y1[i1]  y1[i1m1] ); + ifascend = ( ifxsort && x1[i1] > x1[i1m1] )  + ( !ifxsort && y1[i1] > y1[i1m1] ); + + /* Determine all crossings between the line between the (i11)th + * and i1th vertices of polygon 1 and all edges of polygon 2. + * Retain the lowest and (if ncross ==2) next lowest crossings in + * order of x (or y if ifxsort is false) along the line from i11 + * to i1. */ + + i1m1 = i1  1; + if ( i1m1 < 0 ) + i1m1 += n1; + i2m1 = n2  1; + for ( i2 = 0; i2 < n2; i2++ ) + { + if ( !( x2[i2m1] == x2[i2] && y2[i2m1] == y2[i2] )) + { + ifnotcrossed = notcrossed( &xintersect, &yintersect, + x1[i1m1], y1[i1m1], x1[i1], y1[i1], + x2[i2m1], y2[i2m1], x2[i2], y2[i2] ); + + if ( !ifnotcrossed ) + { + count_crossings++; + if ( count_crossings == 1 ) + { + xcross[0] = xintersect; + ycross[0] = yintersect; + i2cross[0] = i2; + status = 1; + } + else + { + if (( ifxsort && ifascend && xintersect < xcross[0] )  + ( !ifxsort && ifascend && yintersect < ycross[0] )  + ( ifxsort && !ifascend && xintersect >= xcross[0] )  + ( !ifxsort && !ifascend && yintersect >= ycross[0] )) + { + if ( ncross == 2 ) + { + xcross[1] = xcross[0]; + ycross[1] = ycross[0]; + i2cross[1] = i2cross[0]; + status = 2; + } + xcross[0] = xintersect; + ycross[0] = yintersect; + i2cross[0] = i2; + } + else if ( ncross == 2 && ( count_crossings == 2  ( + ( ifxsort && ifascend && xintersect < xcross[1] )  + ( !ifxsort && ifascend && yintersect < ycross[1] )  + ( ifxsort && !ifascend && xintersect >= xcross[1] )  + ( !ifxsort && !ifascend && yintersect >= ycross[1] )))) + { + xcross[1] = xintersect; + ycross[1] = yintersect; + i2cross[1] = i2; + status = 2; + } + } + } + } + i2m1 = i2; + } + return status; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20091228 17:41:52

Revision: 10731 http://plplot.svn.sourceforge.net/plplot/?rev=10731&view=rev Author: airwin Date: 20091228 17:41:45 +0000 (Mon, 28 Dec 2009) Log Message:  Fix cases where split 1/split 2 encompasses all/none of the polygon 2 vertices. This DUSE_FILL_INTERSECTION_POLYGON=ON change finally gives good results for the first four pages of example 25, if the clip limits are shifted around to avoid the case (not implemented yet) where polygon crossings occur at vertices. So it appears I have the basic recursive algorithm working correctly, but there are some details (e.g., polygon crossings at vertices) still to be implemented. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20091227 19:27:12 UTC (rev 10730) +++ trunk/src/plfill.c 20091228 17:41:45 UTC (rev 10731) @@ 1391,7 +1391,8 @@ PLINT i1, i1m1, i1start_new, i2, i2m1, kk, kkstart1, kkstart21, kkstart22,  k, kstart, range1, range21, range22, ncrossed, ncrossed_change, + k, kstart, range1, + range21, range22, ncrossed, ncrossed_change, nsplit1, nsplit2, nsplit2m1; PLINT xintersect[2], yintersect[2], i1intersect[2], i2intersect[2], ifnotcrossed; @@ 1520,7 +1521,6 @@ i1m1 += n1; for ( i1 = i1start; i1 < n1; i1++ ) {  i2m1 = n2  1; ncrossed_change = number_crossings( &xintersect[ncrossed], &yintersect[ncrossed], &i2intersect[ncrossed], 2  ncrossed, @@ 1537,95 +1537,6 @@ /* Have discovered the first two crossings for * polygon 1 at i1 = i1start or above. */  /* Calculate polygon 2 index range in split 1 (the  * split that proceeds beyond the second intersect with  * ascending i2 values). */  range21 = i2intersect[0]  i2intersect[1];  if ( range21 < 0 )  range21 += n2;  /* i2 intersect values range between 0 and n2  1 so  * the the smallest untransformed range21 value is n2  + 1, and the largest untransformed range21 value is  * n2  1. This means the smallest transformed range21  * value is +1. Also, for the untransformed range21  * value of 0, go the "long way" around with  * transformed range21 = n2 so that split 1 is  * guaranteed to have the correct fill_status == +1. */   /* N.B. always avoid special code for range21 == 0  * case because it is no longer necessary since there  * is a guarantee that the first vertex of polygon 1  * (starting with n1  1) that is not near the border  * of polygon 2 will be outside polygon 2. */  if ( 0 && range21 == 0 )  {  /* For this special case the above ascii art  * does not apply, and we must use the  * following diagram instead:  *  *  1 1 ...  *  *  *  * 2???2 X X 2???2  *  * 1 1  *  *  * N.B. no valid split of polygon 2 is  * possible with this diagram. However,  * this diagrem can be classified the same  * as the first diagram, but with the roles  * of polygon 1 and polygon 2 swapped so  * that polygon 1 is the one that gets split  * instead of polygon 2. In fact, swapping  * those two polygons (and starting the new  * polygon 1 just at i2intersect[1] of the old  * polygon 2 to insure the present two  * intersections are immediately (at least for the  * first new polygon segment 1) found and dealt  * with to eliminate the chance of an infinite  * "swap" loop) is the only clean way to deal with  * this situation. */  if (( xsplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ysplit2 = (PLINT *) malloc( n2 * sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  if (( ifsplit1 = (PLINT *) calloc( n1, sizeof ( PLINT ))) == NULL )  {  plexit( "fill_intersection_polygon: Insufficient memory" );  }  kk = i2intersect[1];  for ( i2 = 0; i2 < n2; i2++ )  {  xsplit2[kk] = x2[i2];  ysplit2[kk++] = y2[i2];  if ( kk == n2 )  kk = 0;  }  /* N.B. the positive orientation of split2  * is preserved since the index order is  * the same as that of polygon 2, and by  * assumption that polygon and polygon 1  * have identical positive  * orientations. fill_status does not matter  * here because a split is guaranteed at the  * next recursion level. For definiteness set  * it to zero. */  fill_intersection_polygon(  recursion_depth + 1, ifextrapolygon,  0, fill,  xsplit2, ysplit2, 0, n2,  x1, y1, ifsplit1, n1 );  free( xsplit2 );  free( ysplit2 );  free( ifsplit1 );  return;  } /* New i1start is the same as the current i1 (just * in case there are more crossings to find between * i1m1 and i1.) */ @@ 1645,17 +1556,57 @@ * kkstart1 of polygon 1, and the second intersect. */ kkstart1 = i1intersect[0]; + /* Calculate polygon 2 index range in split 1 (the + * split that proceeds beyond the second intersect with + * ascending i2 values). */ + range21 = i2intersect[0]  i2intersect[1]; + if ( range21 < 0 ) + range21 += n2; + /* i2 intersect values range between 0 and n2  1 so + * the smallest untransformed range21 value is n2 + 1, + * and the largest untransformed range21 value is n2  1. + * This means the smallest transformed range21 value is 0 + * (which occurs only ifi2intersect[0] = i2intersect[1], + * see more commentary for that special case below) while + * the largest transformed range21 value is n2  1. */ + + if ( range21 == 0 ) + { + int ifxsort, ifascend; + /* For this case, the two crossings occur within the same + * polygon 2 boundary segment and if those two crossings + * are in ascending/descending order in i2, then split 1 + * (the split with the positive fill_status) must include + * all/none of the points in polygon 2. */ + i2 = i2intersect[1]; + i2m1 = i2  1; + if ( i2m1 < 0 ) + i2m1 += n2; + + ifxsort = abs( x2[i2]  x2[i2m1] ) > abs( y2[i2]  y2[i2m1] ); + ifascend = ( ifxsort && x2[i2] > x2[i2m1] )  + ( !ifxsort && y2[i2] > y2[i2m1] ); + if (( ifxsort && ifascend && xintersect[0] < xintersect[1] )  + ( !ifxsort && ifascend && yintersect[0] < yintersect[1] )  + ( ifxsort && !ifascend && xintersect[0] >= xintersect[1] )  + ( !ifxsort && !ifascend && yintersect[0] >= yintersect[1] )) + { + range21 = n2; + } + } + kkstart21 = i2intersect[1]; + nsplit1 = 2 + range1 + range21; /* Split 2 of polygon 2 consists of the * boundary + range22 (= n2  range21) points * between kkstart22 (= i2intersect[1]1) and i2intersect[0] in * descending order of polygon 2 indices. */  range22 = n2  range21; + range22 = n2  range21; + /* Starting i2 index of split 2. */ kkstart22 = i2intersect[1]  1; if ( kkstart22 < 0 ) kkstart22 += n2;  nsplit1 = 2 + range1 + range21; nsplit2 = 2 + range1 + range22; if (( xsplit1 = (PLINT *) malloc( nsplit1 * sizeof ( PLINT ))) == NULL ) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20100314 17:18:04

Revision: 10869 http://plplot.svn.sourceforge.net/plplot/?rev=10869&view=rev Author: airwin Date: 20100314 17:17:57 +0000 (Sun, 14 Mar 2010) Log Message:  Remove some of my local work that is not ready for prime time and which was inadvertently included in the commit of David's "Support arbitrary storage of 2D user data" patch. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20100314 17:04:57 UTC (rev 10868) +++ trunk/src/plfill.c 20100314 17:17:57 UTC (rev 10869) @@ 1895,10 +1895,9 @@ status = status  PL_NEAR_PARALLEL; else status = status  PL_PARALLEL;  /* Choice of intersect is arbitrary in this case. Choose A1  * A2, B1, or B2 (in that order) if any of them lie inside or  * near the other line segment. Otherwise, choose the average  * point. */ + /* Choice of intersect is arbitrary in this case. Choose A1, A2, + * B1, or B2 (in that order) if any of them lie inside or near + * the other line segment. Otherwise, choose the average point. */ if ( ( BETW_NBCC( xA1, xB1, xB2 ) && BETW_NBCC( yA1, yB1, yB2 ) ) ) { fxintersect = xA1; @@ 1946,11 +1945,11 @@ * Find out which. */ if ( fabs( fxintersect  xA1 ) <= PL_NBCC && fabs( fyintersect  yA1 ) <= PL_NBCC ) status = status  PL_NEAR_A1;  if ( fabs( fxintersect  xA2 ) <= PL_NBCC && fabs( fyintersect  yA2 ) <= PL_NBCC ) + else if ( fabs( fxintersect  xA2 ) <= PL_NBCC && fabs( fyintersect  yA2 ) <= PL_NBCC ) status = status  PL_NEAR_A2;  if ( fabs( fxintersect  xB1 ) <= PL_NBCC && fabs( fyintersect  yB1 ) <= PL_NBCC ) + else if ( fabs( fxintersect  xB1 ) <= PL_NBCC && fabs( fyintersect  yB1 ) <= PL_NBCC ) status = status  PL_NEAR_B1;  if ( fabs( fxintersect  xB2 ) <= PL_NBCC && fabs( fyintersect  yB2 ) <= PL_NBCC ) + else if ( fabs( fxintersect  xB2 ) <= PL_NBCC && fabs( fyintersect  yB2 ) <= PL_NBCC ) status = status  PL_NEAR_B2; /* N.B. if none of the above conditions hold then status remains at * zero to signal we have a definite crossing. */ @@ 2022,13 +2021,6 @@ int i1m1, i2, i2m1, ifnotcrossed; int ifxsort, ifascend, count_crossings = 0, status = 0; PLINT xintersect, yintersect;  PLINT ifnear_a1, ifnear_a2, ifnear_a, ifnear_b1, ifnear_b2, ifnear_b;  PLINT crossing_index[4];  /* The following are floating point to avoid integer overflows while  * still having exact arithmetic (except for very large values beyond  * the 53 bits of the significand). */  PLFLT xline, yline, xpoint0, ypoint0, xpoint1, ypoint1,  inprod0, inprod1, product_NBCC; i1m1 = i1  1; if ( i1m1 < 0 ) @@ 2061,124 +2053,7 @@ ifnotcrossed = notcrossed( &xintersect, &yintersect, x1[i1m1], y1[i1m1], x1[i1], y1[i1], x2[i2m1], y2[i2m1], x2[i2], y2[i2] );  if ( ifnotcrossed && !( ifnotcrossed & ( PL_PARALLEL  PL_NEAR_PARALLEL  PL_NOT_CROSSED ) ) )  {  /* intersection near at least one of the vertices, but not  * near parallel and not parallel */  ifnear_a1 = ifnotcrossed & PL_NEAR_A1;  ifnear_a2 = ifnotcrossed & PL_NEAR_A2;  ifnear_a = ifnotcrossed & ( PL_NEAR_A1  PL_NEAR_A2 );  ifnear_b1 = ifnotcrossed & PL_NEAR_B1;  ifnear_b2 = ifnotcrossed & PL_NEAR_B2;  ifnear_b = ifnotcrossed & ( PL_NEAR_B1  PL_NEAR_B2 );  if ( ( ifnear_a1 && ifnear_a2 )  ( ifnear_b1 && ifnear_b2 )   !( ifnear_a  ifnear_b ) )  plwarn( "number_crossings; internal logic error" );  if ( ifnear_a1 )  {  crossing_index[0] = i1;  crossing_index[2] = i1m1  1;  if ( crossing_index[2] < 0 )  crossing_index[2] += n1;  }  else if ( ifnear_a2 )  {  crossing_index[0] = i1m1;  crossing_index[2] = i1 + 1;  if ( crossing_index[2] >= n1 )  crossing_index[2] = n1;  }  if ( ifnear_b1 )  {  crossing_index[1] = i2;  crossing_index[3] = i2m1  1;  if ( crossing_index[3] < 0 )  crossing_index[3] += n2;  }  else if ( ifnear_b2 )  {  crossing_index[1] = i2m1;  crossing_index[3] = i2 + 1;  if ( crossing_index[3] >= n2 )  crossing_index[3] = n2;  }  /* Never count ifnear_a1 since that should be a duplicate  * of ifnear_a2 for the next i1 index. */  if ( ifnear_a2 && !ifnear_b )  {  /* Are the points defined by crossing_index[0] (point  * 0) and crossing_index[2] (point 1) on opposite  * sides of the line defined by the polygon 2 edge  * from i2m1 to i2? If so, indicate definite  * crossing. */  /* Use point i2m1 as origin of coordinate system and take  * dot products of both point 0 and point 1 vectors with  * perpendicular to point i2 vector. */  xline = ( y2[i2]  y2[i2m1] );  yline = x2[i2]  x2[i2m1];  xpoint0 = x1[crossing_index[0]]  x2[i2m1];  ypoint0 = y1[crossing_index[0]]  y2[i2m1];  xpoint1 = x1[crossing_index[2]]  x2[i2m1];  ypoint1 = y1[crossing_index[2]]  y2[i2m1];  inprod0 = xpoint0 * xline + ypoint0 * yline;  inprod1 = xpoint1 * xline + ypoint1 * yline;  /* maximum error of product of inner products for  * an error of PL_NBCC of the right sign for  * each original position. */  product_NBCC = 2. * PL_NBCC * (  ( fabs( xpoint0 ) + fabs( xline ) +  fabs( ypoint0 ) + fabs( yline ) ) * fabs( inprod1 ) +  ( fabs( xpoint1 ) + fabs( xline ) +  fabs( ypoint1 ) + fabs( yline ) ) * fabs( inprod0 ) );   /* Not crossed if sine 0 proportional to inprod0 has  * same sign (subject to error) as sine 1 proportional  * to inprod1. */  ifnotcrossed = inprod0 * inprod1 >= product_NBCC;  }  /* Never count ifnear_b1 since that should be a duplicate  * of ifnear_b2 for the next i2 index. */  else if ( ifnear_b2 && !ifnear_a )  {  /* Are the points defined by crossing_index[1] (point  * 0) and crossing_index[3] (point 1) on opposite  * sides of the line defined by the polygon 1 edge  * from i1m1 to i1? If so, indicate definite  * crossing. */   /* Use point i1m1 as origin of coordinate system and take  * dot products of both point 0 and point 1 vectors with  * perpendicular to point i1 vector. */  xline = ( y1[i1]  y1[i1m1] );  yline = x1[i1]  x1[i1m1];  xpoint0 = x2[crossing_index[1]]  x1[i1m1];  ypoint0 = y2[crossing_index[1]]  y1[i1m1];  xpoint1 = x2[crossing_index[3]]  x1[i1m1];  ypoint1 = y2[crossing_index[3]]  y1[i1m1];  inprod0 = xpoint0 * xline + ypoint0 * yline;  inprod1 = xpoint1 * xline + ypoint1 * yline;  /* maximum error of product of inner products for  * an error of PL_NBCC of the right sign for  * each original position. */  product_NBCC = 2. * PL_NBCC * (  ( fabs( xpoint0 ) + fabs( xline ) +  fabs( ypoint0 ) + fabs( yline ) ) * fabs( inprod1 ) +  ( fabs( xpoint1 ) + fabs( xline ) +  fabs( ypoint1 ) + fabs( yline ) ) * fabs( inprod0 ) );   /* Not crossed if sine 0 proportional to inprod0 has  * same sign (subject to error) as sine 1 proportional  * to inprod1. */  ifnotcrossed = inprod0 * inprod1 >= product_NBCC;  }  else if ( ifnear_a && ifnear_b )  {  plwarn( "notcrossing: coincident vertices not implemented yet." );  ifnotcrossed = 0;  }  }  if ( !ifnotcrossed ) { count_crossings++; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <airwin@us...>  20101011 22:53:11

Revision: 11258 http://plplot.svn.sourceforge.net/plplot/?rev=11258&view=rev Author: airwin Date: 20101011 22:53:05 +0000 (Mon, 11 Oct 2010) Log Message:  Fixup minor comment indentation issues. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20101011 22:34:18 UTC (rev 11257) +++ trunk/src/plfill.c 20101011 22:53:05 UTC (rev 11258) @@ 889,11 +889,11 @@ // The cases where the fill region is divided 3/1  // LL LR UR UL // ++ ++ ++ ++ //      \ /  //         // \   /     // ++ ++ ++ ++ + // ++ ++ ++ ++ + //      \ /  + //         + // \   /     + // ++ ++ ++ ++ // // Note when we go the long way around, if the direction is reversed the // three vertices must be visited in the opposite order. @@ 1755,7 +1755,7 @@ // This end phase is reached only if no crossings are found. // If a fill_status of +/ 1 is known, use that to fill or not since // +1 corresponds to all of polygon 2 inside polygon 1 and 1 + // +1 corresponds to all of polygon 2 inside polygon 1 and 1 // corresponds to none of polygon 2 inside polygon 1. if ( fill_status == 1 ) return; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <andrewross@us...>  20110803 14:40:31

Revision: 11846 http://plplot.svn.sourceforge.net/plplot/?rev=11846&view=rev Author: andrewross Date: 20110803 14:40:25 +0000 (Wed, 03 Aug 2011) Log Message:  Slightly cleaner fix to plfill. Also fix a typo bug in plfill3 which would have caused a crash for large numbers of points. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20110803 14:20:40 UTC (rev 11845) +++ trunk/src/plfill.c 20110803 14:40:25 UTC (rev 11846) @@ 134,7 +134,7 @@ { PLINT _xpoly[PL_MAXPOLY], _ypoly[PL_MAXPOLY]; PLINT *xpoly, *ypoly;  PLINT i; + PLINT i, npts; PLFLT xt, yt; if ( plsc>level < 3 ) @@ 147,6 +147,7 @@ plabort( "plfill: Not enough points in object" ); return; } + npts = n; if ( n > PL_MAXPOLY  1 ) { xpoly = (PLINT *) malloc( ( n + 1 ) * sizeof ( PLINT ) ); @@ 181,7 +182,7 @@ plP_plfclp( xpoly, ypoly, n, plsc>clpxmi, plsc>clpxma, plsc>clpymi, plsc>clpyma, plP_fill );  if ( xpoly != _xpoly ) + if ( npts > PL_MAXPOLY  1 ) { free( xpoly ); free( ypoly ); @@ 227,7 +228,7 @@ { tx = (PLFLT *) malloc( ( n + 1 ) * sizeof ( PLFLT ) ); ty = (PLFLT *) malloc( ( n + 1 ) * sizeof ( PLFLT ) );  ty = (PLFLT *) malloc( ( n + 1 ) * sizeof ( PLFLT ) ); + tz = (PLFLT *) malloc( ( n + 1 ) * sizeof ( PLFLT ) ); xpoly = (PLINT *) malloc( ( n + 1 ) * sizeof ( PLINT ) ); ypoly = (PLINT *) malloc( ( n + 1 ) * sizeof ( PLINT ) ); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 
From: <andrewross@us...>  20110915 14:55:10

Revision: 11927 http://plplot.svn.sourceforge.net/plplot/?rev=11927&view=rev Author: andrewross Date: 20110915 14:55:03 +0000 (Thu, 15 Sep 2011) Log Message:  When checking if start / end points of polygon to fill are the same, then check the coordinates after translating into PLplot internal integer format. Prevents spurious postscript file differences due to floating point rounding errors in example 27. Modified Paths:  trunk/src/plfill.c Modified: trunk/src/plfill.c ===================================================================  trunk/src/plfill.c 20110915 08:54:18 UTC (rev 11926) +++ trunk/src/plfill.c 20110915 14:55:03 UTC (rev 11927) @@ 171,12 +171,11 @@ ypoly[i] = plP_wcpcy( yt ); }  if ( x[0] != x[n  1]  y[0] != y[n  1] ) + if ( xpoly[0] != xpoly[n  1]  ypoly[0] != ypoly[n  1] ) { n++;  TRANSFORM( x[0], y[0], &xt, &yt );  xpoly[n  1] = plP_wcpcx( xt );  ypoly[n  1] = plP_wcpcy( yt ); + xpoly[n  1] = xpoly[0]; + ypoly[n  1] = ypoly[0]; } plP_plfclp( xpoly, ypoly, n, plsc>clpxmi, plsc>clpxma, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 