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