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