PLplot  5.11.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
plfill.c
Go to the documentation of this file.
1 // Polygon pattern fill.
2 //
3 // Copyright (C) 2004-2015 Alan W. Irwin
4 // Copyright (C) 2005, 2006, 2007, 2008, 2009 Arjen Markus
5 //
6 // This file is part of PLplot.
7 //
8 // PLplot is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU Library General Public License as published
10 // by the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // PLplot is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU Library General Public License for more details.
17 //
18 // You should have received a copy of the GNU Library General Public License
19 // along with PLplot; if not, write to the Free Software
20 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 //
22 //
23 
24 #include "plplotP.h"
25 
26 #define INSIDE( ix, iy ) ( BETW( ix, xmin, xmax ) && BETW( iy, ymin, ymax ) )
27 
28 #define DTOR ( PI / 180. )
29 #define BINC 50
30 // Near-border comparison criterion (NBCC).
31 #define PL_NBCC 2
32 // Variant of BETW that returns true if between or within PL_NBCC of it.
33 #define BETW_NBCC( ix, ia, ib ) ( ( ( ix ) <= ( ia + PL_NBCC ) && ( ix ) >= ( ib - PL_NBCC ) ) || ( ( ix ) >= ( ia - PL_NBCC ) && ( ix ) <= ( ib + PL_NBCC ) ) )
34 
35 // Status codes ORed together in the return from notcrossed.
36 // PL_NOT_CROSSED is set whenever the two line segments definitely
37 // (i.e., intersection not near the ends or completely apart)
38 // do not cross each other.
39 //
40 // PL_NEAR_A1 is set if the intersect is near (+/- PL_NBCC) the first
41 // end of line segment A.
42 //
43 // PL_NEAR_A2 is set if the intersect is near (+/- PL_NBCC) the second
44 // end of line segment A.
45 //
46 // PL_NEAR_B1 is set if the intersect is near (+/- PL_NBCC) the first
47 // end of line segment B.
48 //
49 // PL_NEAR_B2 is set if the intersect is near (+/- PL_NBCC) the second
50 // end of line segment B.
51 //
52 // PL_NEAR_PARALLEL is set if the two line segments are nearly parallel
53 // with each other, i.e., a change in one of the coordinates of up to
54 // (+/- PL_NBCC) would render them exactly parallel.
55 //
56 // PL_PARALLEL is set if the two line segments are exactly parallel
57 // with each other.
58 //
60 {
62  PL_NEAR_A1 = 0x2,
63  PL_NEAR_A2 = 0x4,
64  PL_NEAR_B1 = 0x8,
65  PL_NEAR_B2 = 0x10,
67  PL_PARALLEL = 0x40
68 };
69 
70 struct point
71 {
72  PLINT x, y;
73 };
75 
76 // Static function prototypes
77 
78 static int
79 compar( const void *, const void * );
80 
81 static void
82 addcoord( PLINT, PLINT );
83 
84 static void
85 tran( PLINT *, PLINT *, PLFLT, PLFLT );
86 
87 static void
89 
90 static int
91 notpointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp );
92 
93 static int
94 circulation( PLINT *x, PLINT *y, PLINT npts );
95 
96 #ifdef USE_FILL_INTERSECTION_POLYGON
97 static void
98 fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon,
99  PLINT fill_status,
100  void ( *fill )( short *, short *, PLINT ),
101  const PLINT *x1, const PLINT *y1,
102  PLINT i1start, PLINT n1,
103  const PLINT *x2, const PLINT *y2,
104  const PLINT *if2, PLINT n2 );
105 
106 static int
107 positive_orientation( PLINT n, const PLINT *x, const PLINT *y );
108 
109 static int
110 number_crossings( PLINT *xcross, PLINT *ycross, PLINT *i2cross, PLINT ncross,
111  PLINT i1, PLINT n1, const PLINT *x1, const PLINT *y1,
112  PLINT n2, const PLINT *x2, const PLINT *y2 );
113 #endif
114 
115 static int
116 notcrossed( PLINT *xintersect, PLINT *yintersect,
117  PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2,
118  PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 );
119 
120 //--------------------------------------------------------------------------
121 // void plfill()
122 //
123 // Pattern fills the polygon bounded by the input points.
124 // For a number of vertices greater than PL_MAXPOLY-1, memory is managed via
125 // malloc/free. Otherwise statically defined arrays of length PL_MAXPOLY
126 // are used.
127 // The final point is explicitly added if it doesn't match up to the first,
128 // to prevent clipping problems.
129 //--------------------------------------------------------------------------
130 
131 void
132 c_plfill( PLINT n, const PLFLT *x, const PLFLT *y )
133 {
134  PLINT _xpoly[PL_MAXPOLY], _ypoly[PL_MAXPOLY];
135  PLINT *xpoly, *ypoly;
136  PLINT i, npts;
137  PLFLT xt, yt;
138 
139  if ( plsc->level < 3 )
140  {
141  plabort( "plfill: Please set up window first" );
142  return;
143  }
144  if ( n < 3 )
145  {
146  plabort( "plfill: Not enough points in object" );
147  return;
148  }
149  npts = n;
150  if ( n > PL_MAXPOLY - 1 )
151  {
152  xpoly = (PLINT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLINT ) );
153  ypoly = (PLINT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLINT ) );
154 
155  if ( ( xpoly == NULL ) || ( ypoly == NULL ) )
156  {
157  plexit( "plfill: Insufficient memory for large polygon" );
158  }
159  }
160  else
161  {
162  xpoly = _xpoly;
163  ypoly = _ypoly;
164  }
165 
166  for ( i = 0; i < n; i++ )
167  {
168  TRANSFORM( x[i], y[i], &xt, &yt );
169  xpoly[i] = plP_wcpcx( xt );
170  ypoly[i] = plP_wcpcy( yt );
171  }
172 
173  if ( xpoly[0] != xpoly[n - 1] || ypoly[0] != ypoly[n - 1] )
174  {
175  n++;
176  xpoly[n - 1] = xpoly[0];
177  ypoly[n - 1] = ypoly[0];
178  }
179 
180  plP_plfclp( xpoly, ypoly, n, plsc->clpxmi, plsc->clpxma,
181  plsc->clpymi, plsc->clpyma, plP_fill );
182 
183  if ( npts > PL_MAXPOLY - 1 )
184  {
185  free( xpoly );
186  free( ypoly );
187  }
188 }
189 
190 //--------------------------------------------------------------------------
191 // void plfill3()
192 //
193 // Pattern fills the polygon in 3d bounded by the input points.
194 // For a number of vertices greater than PL_MAXPOLY-1, memory is managed via
195 // malloc/free. Otherwise statically defined arrays of length PL_MAXPOLY
196 // are used.
197 // The final point is explicitly added if it doesn't match up to the first,
198 // to prevent clipping problems.
199 //--------------------------------------------------------------------------
200 
201 void
202 c_plfill3( PLINT n, const PLFLT *x, const PLFLT *y, const PLFLT *z )
203 {
204  PLFLT _tx[PL_MAXPOLY], _ty[PL_MAXPOLY], _tz[PL_MAXPOLY];
205  PLFLT *tx, *ty, *tz;
206  PLFLT *V[3];
207  PLINT _xpoly[PL_MAXPOLY], _ypoly[PL_MAXPOLY];
208  PLINT *xpoly, *ypoly;
209  PLINT i;
210  PLINT npts;
211  PLFLT xmin, xmax, ymin, ymax, zmin, zmax, zscale;
212 
213  if ( plsc->level < 3 )
214  {
215  plabort( "plfill3: Please set up window first" );
216  return;
217  }
218  if ( n < 3 )
219  {
220  plabort( "plfill3: Not enough points in object" );
221  return;
222  }
223 
224  npts = n;
225  if ( n > PL_MAXPOLY - 1 )
226  {
227  tx = (PLFLT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLFLT ) );
228  ty = (PLFLT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLFLT ) );
229  tz = (PLFLT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLFLT ) );
230  xpoly = (PLINT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLINT ) );
231  ypoly = (PLINT *) malloc( (size_t) ( n + 1 ) * sizeof ( PLINT ) );
232 
233  if ( ( tx == NULL ) || ( ty == NULL ) || ( tz == NULL ) ||
234  ( xpoly == NULL ) || ( ypoly == NULL ) )
235  {
236  plexit( "plfill3: Insufficient memory for large polygon" );
237  }
238  }
239  else
240  {
241  tx = _tx;
242  ty = _ty;
243  tz = _tz;
244  xpoly = _xpoly;
245  ypoly = _ypoly;
246  }
247 
248  plP_gdom( &xmin, &xmax, &ymin, &ymax );
249  plP_grange( &zscale, &zmin, &zmax );
250 
251  // copy the vertices so we can clip without corrupting the input
252  for ( i = 0; i < n; i++ )
253  {
254  tx[i] = x[i]; ty[i] = y[i]; tz[i] = z[i];
255  }
256  if ( tx[0] != tx[n - 1] || ty[0] != ty[n - 1] || tz[0] != tz[n - 1] )
257  {
258  n++;
259  tx[n - 1] = tx[0]; ty[n - 1] = ty[0]; tz[n - 1] = tz[0];
260  }
261  V[0] = tx; V[1] = ty; V[2] = tz;
262  n = plP_clip_poly( n, V, 0, 1, -xmin );
263  n = plP_clip_poly( n, V, 0, -1, xmax );
264  n = plP_clip_poly( n, V, 1, 1, -ymin );
265  n = plP_clip_poly( n, V, 1, -1, ymax );
266  n = plP_clip_poly( n, V, 2, 1, -zmin );
267  n = plP_clip_poly( n, V, 2, -1, zmax );
268  for ( i = 0; i < n; i++ )
269  {
270  xpoly[i] = plP_wcpcx( plP_w3wcx( tx[i], ty[i], tz[i] ) );
271  ypoly[i] = plP_wcpcy( plP_w3wcy( tx[i], ty[i], tz[i] ) );
272  }
273 
274 // AWI: in the past we have used
275 // plP_fill(xpoly, ypoly, n);
276 // here, but our educated guess is this fill should be done via the clipping
277 // interface instead as below.
278 // No example tests this code so one of our users will end up inadvertently
279 // testing this for us.
280 //
281 // jc: I have checked, and both versions does give the same result, i.e., clipping
282 // to the window boundaries. The reason is that the above plP_clip_poly() does
283 // the clipping. To check this, is enough to diminish the x/y/z min/max arguments in
284 // plw3d() in x08c. But let's keep it, although 10% slower...
285 //
286  plP_plfclp( xpoly, ypoly, n, plsc->clpxmi, plsc->clpxma,
287  plsc->clpymi, plsc->clpyma, plP_fill );
288 
289 // If the original number of points is large, then free the arrays
290  if ( npts > PL_MAXPOLY - 1 )
291  {
292  free( tx );
293  free( ty );
294  free( tz );
295  free( xpoly );
296  free( ypoly );
297  }
298 }
299 
300 //--------------------------------------------------------------------------
301 // void plfill_soft()
302 //
303 // Pattern fills in software the polygon bounded by the input points.
304 //--------------------------------------------------------------------------
305 
306 void
307 plfill_soft( short *x, short *y, PLINT n )
308 {
309  PLINT i, j;
310  PLINT xp1, yp1, xp2, yp2, xp3, yp3;
311  PLINT k, dinc;
312  PLFLT ci, si;
314  double temp;
315 
316  buffersize = 2 * BINC;
317  buffer = (PLINT *) malloc( (size_t) buffersize * sizeof ( PLINT ) );
318  if ( !buffer )
319  {
320  plabort( "plfill: Out of memory" );
321  return;
322  }
323 
324  //do not write the hatching lines to the buffer as we have already
325  //written the fill to the buffer
326  plbuf_write = plsc->plbuf_write;
327  plsc->plbuf_write = FALSE;
328 // Loop over sets of lines in pattern
329 
330  for ( k = 0; k < plsc->nps; k++ )
331  {
332  bufferleng = 0;
333 
334  temp = DTOR * plsc->inclin[k] * 0.1;
335  si = sin( temp ) * plsc->ypmm;
336  ci = cos( temp ) * plsc->xpmm;
337 
338  // normalize: 1 = si*si + ci*ci
339 
340  temp = sqrt( (double) ( si * si + ci * ci ) );
341  si /= temp;
342  ci /= temp;
343 
344  dinc = (PLINT) ( plsc->delta[k] * SSQR( plsc->ypmm * ABS( ci ),
345  plsc->xpmm * ABS( si ) ) / 1000. );
346 
347  if ( dinc < 0 )
348  dinc = -dinc;
349  if ( dinc == 0 )
350  dinc = 1;
351 
352  xp1 = x[n - 2];
353  yp1 = y[n - 2];
354  tran( &xp1, &yp1, (PLFLT) ci, (PLFLT) si );
355 
356  xp2 = x[n - 1];
357  yp2 = y[n - 1];
358  tran( &xp2, &yp2, (PLFLT) ci, (PLFLT) si );
359 
360 // Loop over points in polygon
361 
362  for ( i = 0; i < n; i++ )
363  {
364  xp3 = x[i];
365  yp3 = y[i];
366  tran( &xp3, &yp3, (PLFLT) ci, (PLFLT) si );
367  buildlist( xp1, yp1, xp2, yp2, xp3, yp3, dinc );
368  xp1 = xp2;
369  yp1 = yp2;
370  xp2 = xp3;
371  yp2 = yp3;
372  }
373 
374 // Sort list by y then x
375 
376  qsort( (void *) buffer, (size_t) bufferleng / 2,
377  (size_t) sizeof ( struct point ), compar );
378 
379 // OK, now do the hatching
380 
381  i = 0;
382 
383  while ( i < bufferleng )
384  {
385  xp1 = buffer[i];
386  yp1 = buffer[i + 1];
387  i += 2;
388  xp2 = xp1;
389  yp2 = yp1;
390  tran( &xp1, &yp1, (PLFLT) ci, (PLFLT) ( -si ) );
391  plP_movphy( xp1, yp1 );
392  xp1 = buffer[i];
393  yp1 = buffer[i + 1];
394  i += 2;
395  if ( yp2 != yp1 )
396  {
397  fprintf( stderr, "plfill: oh oh we are lost\n" );
398  for ( j = 0; j < bufferleng; j += 2 )
399  {
400  fprintf( stderr, "plfill: %d %d\n",
401  (int) buffer[j], (int) buffer[j + 1] );
402  }
403  continue; // Uh oh we're lost
404  }
405  tran( &xp1, &yp1, (PLFLT) ci, (PLFLT) ( -si ) );
406  plP_draphy( xp1, yp1 );
407  }
408  }
409  //reinstate the buffer writing parameter and free memory
410  plsc->plbuf_write = plbuf_write;
411  free( (void *) buffer );
412 }
413 
414 //--------------------------------------------------------------------------
415 // Utility functions
416 //--------------------------------------------------------------------------
417 
418 void
419 tran( PLINT *a, PLINT *b, PLFLT c, PLFLT d )
420 {
421  PLINT ta, tb;
422 
423  ta = *a;
424  tb = *b;
425 
426  *a = (PLINT) floor( (double) ( ta * c + tb * d + 0.5 ) );
427  *b = (PLINT) floor( (double) ( tb * c - ta * d + 0.5 ) );
428 }
429 
430 void
431 buildlist( PLINT xp1, PLINT yp1, PLINT xp2, PLINT yp2, PLINT PL_UNUSED( xp3 ), PLINT yp3,
432  PLINT dinc )
433 {
434  PLINT min_y, max_y;
435  PLINT dx, dy, cstep, nstep, ploty, plotx;
436 
437  dx = xp2 - xp1;
438  dy = yp2 - yp1;
439 
440  if ( dy == 0 )
441  {
442  if ( yp2 > yp3 && ( ( yp2 % dinc ) == 0 ) )
443  addcoord( xp2, yp2 );
444  return;
445  }
446 
447  if ( dy > 0 )
448  {
449  cstep = 1;
450  min_y = yp1;
451  max_y = yp2;
452  }
453  else
454  {
455  cstep = -1;
456  min_y = yp2;
457  max_y = yp1;
458  }
459 
460  nstep = ( yp3 > yp2 ? 1 : -1 );
461  if ( yp3 == yp2 )
462  nstep = 0;
463 
464  // Build coordinate list
465 
466  ploty = ( min_y / dinc ) * dinc;
467  if ( ploty < min_y )
468  ploty += dinc;
469 
470  for (; ploty <= max_y; ploty += dinc )
471  {
472  if ( ploty == yp1 )
473  continue;
474  if ( ploty == yp2 )
475  {
476  if ( cstep == -nstep )
477  continue;
478  if ( yp2 == yp3 && yp1 > yp2 )
479  continue;
480  }
481  plotx = xp1 + (PLINT) floor( ( (double) ( ploty - yp1 ) * dx ) / dy + 0.5 );
482  addcoord( plotx, ploty );
483  }
484 }
485 
486 void
487 addcoord( PLINT xp1, PLINT yp1 )
488 {
489  PLINT *temp;
490 
491  if ( bufferleng + 2 > buffersize )
492  {
493  buffersize += 2 * BINC;
494  temp = (PLINT *) realloc( (void *) buffer,
495  (size_t) buffersize * sizeof ( PLINT ) );
496  if ( !temp )
497  {
498  free( (void *) buffer );
499  plexit( "plfill: Out of memory!" );
500  }
501  buffer = temp;
502  }
503 
504  buffer[bufferleng++] = xp1;
505  buffer[bufferleng++] = yp1;
506 }
507 
508 int
509 compar( const void *pnum1, const void *pnum2 )
510 {
511  const struct point *pnt1, *pnt2;
512 
513  pnt1 = (const struct point *) pnum1;
514  pnt2 = (const struct point *) pnum2;
515 
516  if ( pnt1->y < pnt2->y )
517  return -1;
518  else if ( pnt1->y > pnt2->y )
519  return 1;
520 
521  // Only reach here if y coords are equal, so sort by x
522 
523  if ( pnt1->x < pnt2->x )
524  return -1;
525  else if ( pnt1->x > pnt2->x )
526  return 1;
527 
528  return 0;
529 }
530 
531 //--------------------------------------------------------------------------
532 // void plP_plfclp()
533 //
534 // Fills a polygon within the clip limits.
535 //--------------------------------------------------------------------------
536 
537 void
540  void ( *draw )( short *, short *, PLINT ) )
541 {
542 #ifdef USE_FILL_INTERSECTION_POLYGON
543  PLINT *x10, *y10, *x1, *y1, *if1, i1start = 0, i, im1, n1, n1m1,
544  ifnotpointinpolygon;
545  PLINT x2[4] = { xmin, xmax, xmax, xmin };
546  PLINT y2[4] = { ymin, ymin, ymax, ymax };
547  PLINT if2[4] = { 0, 0, 0, 0 };
548  PLINT n2 = 4;
549 
550  // Must have at least 3 points and draw() specified
551  if ( npts < 3 || !draw )
552  return;
553 
554  if ( ( x10 = (PLINT *) malloc( (size_t) npts * sizeof ( PLINT ) ) ) == NULL )
555  {
556  plexit( "plP_plfclp: Insufficient memory" );
557  }
558  if ( ( y10 = (PLINT *) malloc( (size_t) npts * sizeof ( PLINT ) ) ) == NULL )
559  {
560  plexit( "plP_plfclp: Insufficient memory" );
561  }
562  // Polygon 2 obviously has no dups nor two consective segments that
563  // are parallel, but get rid of those type of segments in polygon 1
564  // if they exist.
565 
566  im1 = npts - 1;
567  n1 = 0;
568  for ( i = 0; i < npts; i++ )
569  {
570  if ( !( x[i] == x[im1] && y[i] == y[im1] ) )
571  {
572  x10[n1] = x[i];
573  y10[n1++] = y[i];
574  }
575  im1 = i;
576  }
577 
578  // Must have at least three points that satisfy the above criteria.
579  if ( n1 < 3 )
580  {
581  free( x10 );
582  free( y10 );
583  return;
584  }
585 
586  // Polygon 2 obviously has a positive orientation (i.e., as you
587  // ascend in index along the boundary, the points just adjacent to
588  // the boundary and on the left are interior points for the
589  // polygon), but enforce this condition demanded by
590  // fill_intersection_polygon for polygon 1 as well.
591  if ( positive_orientation( n1, x10, y10 ) )
592  {
593  x1 = x10;
594  y1 = y10;
595  }
596  else
597  {
598  if ( ( x1 = (PLINT *) malloc( (size_t) n1 * sizeof ( PLINT ) ) ) == NULL )
599  {
600  plexit( "plP_plfclp: Insufficient memory" );
601  }
602  if ( ( y1 = (PLINT *) malloc( (size_t) n1 * sizeof ( PLINT ) ) ) == NULL )
603  {
604  plexit( "plP_plfclp: Insufficient memory" );
605  }
606  n1m1 = n1 - 1;
607  for ( i = 0; i < n1; i++ )
608  {
609  x1[n1m1 - i] = x10[i];
610  y1[n1m1 - i] = y10[i];
611  }
612  free( x10 );
613  free( y10 );
614  }
615 
616  // Insure that the first vertex of polygon 1 (starting with n1 -
617  // 1) that is not on the border of polygon 2 is definitely outside
618  // polygon 2.
619  im1 = n1 - 1;
620  for ( i = 0; i < n1; i++ )
621  {
622  if ( ( ifnotpointinpolygon =
623  notpointinpolygon( n2, x2, y2, x1[im1], y1[im1] ) ) != 1 )
624  break;
625  im1 = i;
626  }
627 
628  if ( ifnotpointinpolygon )
629  fill_intersection_polygon( 0, 0, 0, draw, x1, y1, i1start, n1, x2, y2, if2, n2 );
630  else
631  {
632  if ( ( if1 = (PLINT *) calloc( n1, sizeof ( PLINT ) ) ) == NULL )
633  {
634  plexit( "plP_plfclp: Insufficient memory" );
635  }
636  fill_intersection_polygon( 0, 0, 0, draw, x2, y2, i1start, n2, x1, y1, if1, n1 );
637  free( if1 );
638  }
639  free( x1 );
640  free( y1 );
641  return;
642 }
643 #else // USE_FILL_INTERSECTION_POLYGON
644 
645  PLINT i, x1, x2, y1, y2;
646  int iclp = 0, iout = 2;
647  short _xclp[2 * PL_MAXPOLY + 2], _yclp[2 * PL_MAXPOLY + 2];
648  short *xclp = NULL, *yclp = NULL;
649  int drawable;
650  int crossed_xmin1 = 0, crossed_xmax1 = 0;
651  int crossed_ymin1 = 0, crossed_ymax1 = 0;
652  int crossed_xmin2 = 0, crossed_xmax2 = 0;
653  int crossed_ymin2 = 0, crossed_ymax2 = 0;
654  int crossed_up = 0, crossed_down = 0;
655  int crossed_left = 0, crossed_right = 0;
656  int inside_lb;
657  int inside_lu;
658  int inside_rb;
659  int inside_ru;
660 
661  // Must have at least 3 points and draw() specified
662  if ( npts < 3 || !draw )
663  return;
664 
665  if ( npts < PL_MAXPOLY )
666  {
667  xclp = _xclp;
668  yclp = _yclp;
669  }
670  else
671  {
672  if ( ( ( xclp = (short *) malloc( (size_t) ( 2 * npts + 2 ) * sizeof ( short ) ) ) == NULL ) ||
673  ( ( yclp = (short *) malloc( (size_t) ( 2 * npts + 2 ) * sizeof ( short ) ) ) == NULL ) )
674  {
675  plexit( "plP_plfclp: Insufficient memory" );
676  }
677  }
678  inside_lb = !notpointinpolygon( npts, x, y, xmin, ymin );
679  inside_lu = !notpointinpolygon( npts, x, y, xmin, ymax );
680  inside_rb = !notpointinpolygon( npts, x, y, xmax, ymin );
681  inside_ru = !notpointinpolygon( npts, x, y, xmax, ymax );
682 
683  for ( i = 0; i < npts - 1; i++ )
684  {
685  x1 = x[i]; x2 = x[i + 1];
686  y1 = y[i]; y2 = y[i + 1];
687 
688  drawable = ( INSIDE( x1, y1 ) && INSIDE( x2, y2 ) );
689  if ( !drawable )
690  drawable = !plP_clipline( &x1, &y1, &x2, &y2,
691  xmin, xmax, ymin, ymax );
692 
693  if ( drawable )
694  {
695  // Boundary crossing condition -- coming in.
696  crossed_xmin2 = ( x1 == xmin ); crossed_xmax2 = ( x1 == xmax );
697  crossed_ymin2 = ( y1 == ymin ); crossed_ymax2 = ( y1 == ymax );
698 
699  crossed_left = ( crossed_left || crossed_xmin2 );
700  crossed_right = ( crossed_right || crossed_xmax2 );
701  crossed_down = ( crossed_down || crossed_ymin2 );
702  crossed_up = ( crossed_up || crossed_ymax2 );
703  iout = iclp + 2;
704  // If the first segment, just add it.
705 
706  if ( iclp == 0 )
707  {
708  xclp[iclp] = (short) x1; yclp[iclp] = (short) y1; iclp++;
709  xclp[iclp] = (short) x2; yclp[iclp] = (short) y2; iclp++;
710  }
711 
712  // Not first point. If first point of this segment matches up to the
713  // previous point, just add it.
714 
715  else if ( x1 == (int) xclp[iclp - 1] && y1 == (int) yclp[iclp - 1] )
716  {
717  xclp[iclp] = (short) x2; yclp[iclp] = (short) y2; iclp++;
718  }
719 
720  // Otherwise, we need to add both points, to connect the points in the
721  // polygon along the clip boundary. If we encircled a corner, we have
722  // to add that first.
723  //
724 
725  else
726  {
727  // Treat the case where we encircled two corners:
728  // Construct a polygon out of the subset of vertices
729  // Note that the direction is important too when adding
730  // the extra points
731  xclp[iclp + 1] = (short) x2; yclp[iclp + 1] = (short) y2;
732  xclp[iclp + 2] = (short) x1; yclp[iclp + 2] = (short) y1;
733  iout = iout - iclp + 1;
734  // Upper two
735  if ( ( ( crossed_xmin1 && crossed_xmax2 ) ||
736  ( crossed_xmin2 && crossed_xmax1 ) ) &&
737  inside_lu )
738  {
739  if ( crossed_xmin1 )
740  {
741  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
742  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
743  }
744  else
745  {
746  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
747  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
748  }
749  }
750  // Lower two
751  else if ( ( ( crossed_xmin1 && crossed_xmax2 ) ||
752  ( crossed_xmin2 && crossed_xmax1 ) ) &&
753  inside_lb )
754  {
755  if ( crossed_xmin1 )
756  {
757  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
758  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
759  }
760  else
761  {
762  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
763  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
764  }
765  }
766  // Left two
767  else if ( ( ( crossed_ymin1 && crossed_ymax2 ) ||
768  ( crossed_ymin2 && crossed_ymax1 ) ) &&
769  inside_lb )
770  {
771  if ( crossed_ymin1 )
772  {
773  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
774  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
775  }
776  else
777  {
778  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
779  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
780  }
781  }
782  // Right two
783  else if ( ( ( crossed_ymin1 && crossed_ymax2 ) ||
784  ( crossed_ymin2 && crossed_ymax1 ) ) &&
785  inside_rb )
786  {
787  if ( crossed_ymin1 )
788  {
789  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
790  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
791  }
792  else
793  {
794  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
795  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
796  }
797  }
798  // Now the case where we encircled one corner
799  // Lower left
800  else if ( ( crossed_xmin1 && crossed_ymin2 ) ||
801  ( crossed_ymin1 && crossed_xmin2 ) )
802  {
803  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
804  }
805  // Lower right
806  else if ( ( crossed_xmax1 && crossed_ymin2 ) ||
807  ( crossed_ymin1 && crossed_xmax2 ) )
808  {
809  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
810  }
811  // Upper left
812  else if ( ( crossed_xmin1 && crossed_ymax2 ) ||
813  ( crossed_ymax1 && crossed_xmin2 ) )
814  {
815  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
816  }
817  // Upper right
818  else if ( ( crossed_xmax1 && crossed_ymax2 ) ||
819  ( crossed_ymax1 && crossed_xmax2 ) )
820  {
821  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
822  }
823 
824  // Now add current segment.
825  xclp[iclp] = (short) x1; yclp[iclp] = (short) y1; iclp++;
826  xclp[iclp] = (short) x2; yclp[iclp] = (short) y2; iclp++;
827  }
828 
829  // Boundary crossing condition -- going out.
830  crossed_xmin1 = ( x2 == xmin ); crossed_xmax1 = ( x2 == xmax );
831  crossed_ymin1 = ( y2 == ymin ); crossed_ymax1 = ( y2 == ymax );
832  }
833  }
834 
835  // Limit case - all vertices are outside of bounding box. So just fill entire
836  // box, *if* the bounding box is completely encircled.
837  //
838  if ( iclp == 0 )
839  {
840  if ( inside_lb )
841  {
842  xclp[0] = (short) xmin; yclp[0] = (short) ymin;
843  xclp[1] = (short) xmax; yclp[1] = (short) ymin;
844  xclp[2] = (short) xmax; yclp[2] = (short) ymax;
845  xclp[3] = (short) xmin; yclp[3] = (short) ymax;
846  xclp[4] = (short) xmin; yclp[4] = (short) ymin;
847  ( *draw )( xclp, yclp, 5 );
848 
849  if ( xclp != _xclp )
850  {
851  free( xclp );
852  free( yclp );
853  }
854 
855  return;
856  }
857  }
858 
859  // Now handle cases where fill polygon intersects two sides of the box
860 
861  if ( iclp >= 2 )
862  {
863  int debug = 0;
864  int dir = circulation( x, y, npts );
865  if ( debug )
866  {
867  if ( ( xclp[0] == (short) xmin && xclp[iclp - 1] == (short) xmax ) ||
868  ( xclp[0] == (short) xmax && xclp[iclp - 1] == (short) xmin ) ||
869  ( yclp[0] == (short) ymin && yclp[iclp - 1] == (short) ymax ) ||
870  ( yclp[0] == (short) ymax && yclp[iclp - 1] == (short) ymin ) ||
871  ( xclp[0] == (short) xmin && yclp[iclp - 1] == (short) ymin ) ||
872  ( yclp[0] == (short) ymin && xclp[iclp - 1] == (short) xmin ) ||
873  ( xclp[0] == (short) xmax && yclp[iclp - 1] == (short) ymin ) ||
874  ( yclp[0] == (short) ymin && xclp[iclp - 1] == (short) xmax ) ||
875  ( xclp[0] == (short) xmax && yclp[iclp - 1] == (short) ymax ) ||
876  ( yclp[0] == (short) ymax && xclp[iclp - 1] == (short) xmax ) ||
877  ( xclp[0] == (short) xmin && yclp[iclp - 1] == (short) ymax ) ||
878  ( yclp[0] == (short) ymax && xclp[iclp - 1] == (short) xmin ) )
879  {
880  printf( "dir=%d, clipped points:\n", dir );
881  for ( i = 0; i < iclp; i++ )
882  printf( " x[%d]=%hd y[%d]=%hd", i, xclp[i], i, yclp[i] );
883  printf( "\n" );
884  printf( "pre-clipped points:\n" );
885  for ( i = 0; i < npts; i++ )
886  printf( " x[%d]=%d y[%d]=%d", i, x[i], i, y[i] );
887  printf( "\n" );
888  }
889  }
890 
891  // The cases where the fill region is divided 2/2
892  // Divided horizontally
893  if ( xclp[0] == (short) xmin && xclp[iclp - 1] == (short) xmax )
894  {
895  if ( dir > 0 )
896  {
897  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
898  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
899  }
900  else
901  {
902  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
903  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
904  }
905  }
906  else if ( xclp[0] == (short) xmax && xclp[iclp - 1] == (short) xmin )
907  {
908  if ( dir > 0 )
909  {
910  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
911  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
912  }
913  else
914  {
915  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
916  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
917  }
918  }
919 
920  // Divided vertically
921  else if ( yclp[0] == (short) ymin && yclp[iclp - 1] == (short) ymax )
922  {
923  if ( dir > 0 )
924  {
925  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
926  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
927  }
928  else
929  {
930  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
931  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
932  }
933  }
934  else if ( yclp[0] == (short) ymax && yclp[iclp - 1] == (short) ymin )
935  {
936  if ( dir > 0 )
937  {
938  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
939  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
940  }
941  else
942  {
943  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
944  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
945  }
946  }
947 
948  // The cases where the fill region is divided 3/1 --
949  // LL LR UR UL
950  // +-----+ +-----+ +-----+ +-----+
951  // | | | | | \| |/ |
952  // | | | | | | | |
953  // |\ | | /| | | | |
954  // +-----+ +-----+ +-----+ +-----+
955  //
956  // Note when we go the long way around, if the direction is reversed the
957  // three vertices must be visited in the opposite order.
958  //
959  // LL, short way around
960  else if ( ( xclp[0] == (short) xmin && yclp[iclp - 1] == (short) ymin && dir < 0 ) ||
961  ( yclp[0] == (short) ymin && xclp[iclp - 1] == (short) xmin && dir > 0 ) )
962  {
963  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
964  }
965  // LL, long way around, counterclockwise
966  else if ( ( xclp[0] == (short) xmin && yclp[iclp - 1] == (short) ymin && dir > 0 ) )
967  {
968  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
969  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
970  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
971  }
972  // LL, long way around, clockwise
973  else if ( ( yclp[0] == ymin && xclp[iclp - 1] == xmin && dir < 0 ) )
974  {
975  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
976  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
977  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
978  }
979  // LR, short way around
980  else if ( ( xclp[0] == (short) xmax && yclp[iclp - 1] == (short) ymin && dir > 0 ) ||
981  ( yclp[0] == (short) ymin && xclp[iclp - 1] == (short) xmax && dir < 0 ) )
982  {
983  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
984  }
985  // LR, long way around, counterclockwise
986  else if ( yclp[0] == (short) ymin && xclp[iclp - 1] == (short) xmax && dir > 0 )
987  {
988  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
989  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
990  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
991  }
992  // LR, long way around, clockwise
993  else if ( xclp[0] == (short) xmax && yclp[iclp - 1] == (short) ymin && dir < 0 )
994  {
995  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
996  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
997  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
998  }
999  // UR, short way around
1000  else if ( ( xclp[0] == (short) xmax && yclp[iclp - 1] == (short) ymax && dir < 0 ) ||
1001  ( yclp[0] == (short) ymax && xclp[iclp - 1] == (short) xmax && dir > 0 ) )
1002  {
1003  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
1004  }
1005  // UR, long way around, counterclockwise
1006  else if ( xclp[0] == (short) xmax && yclp[iclp - 1] == (short) ymax && dir > 0 )
1007  {
1008  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
1009  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
1010  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
1011  }
1012  // UR, long way around, clockwise
1013  else if ( yclp[0] == (short) ymax && xclp[iclp - 1] == (short) xmax && dir < 0 )
1014  {
1015  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
1016  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
1017  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
1018  }
1019  // UL, short way around
1020  else if ( ( xclp[0] == (short) xmin && yclp[iclp - 1] == (short) ymax && dir > 0 ) ||
1021  ( yclp[0] == (short) ymax && xclp[iclp - 1] == (short) xmin && dir < 0 ) )
1022  {
1023  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymax; iclp++;
1024  }
1025  // UL, long way around, counterclockwise
1026  else if ( yclp[0] == (short) ymax && xclp[iclp - 1] == (short) xmin && dir > 0 )
1027  {
1028  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
1029  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
1030  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
1031  }
1032  // UL, long way around, clockwise
1033  else if ( xclp[0] == (short) xmin && yclp[iclp - 1] == (short) ymax && dir < 0 )
1034  {
1035  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymax; iclp++;
1036  xclp[iclp] = (short) xmax; yclp[iclp] = (short) ymin; iclp++;
1037  xclp[iclp] = (short) xmin; yclp[iclp] = (short) ymin; iclp++;
1038  }
1039  }
1040 
1041  // Check for the case that only one side has been crossed
1042  // (AM) Just checking a single point turns out not to be
1043  // enough, apparently the crossed_*1 and crossed_*2 variables
1044  // are not quite what I expected.
1045  //
1046  if ( inside_lb + inside_rb + inside_lu + inside_ru == 4 )
1047  {
1048  int dir = circulation( x, y, npts );
1049  PLINT xlim[4], ylim[4];
1050  int insert = -99;
1051  int incr = -99;
1052 
1053  xlim[0] = xmin; ylim[0] = ymin;
1054  xlim[1] = xmax; ylim[1] = ymin;
1055  xlim[2] = xmax; ylim[2] = ymax;
1056  xlim[3] = xmin; ylim[3] = ymax;
1057 
1058  if ( crossed_left + crossed_right + crossed_down + crossed_up == 1 )
1059  {
1060  if ( dir > 0 )
1061  {
1062  incr = 1;
1063  insert = 0 * crossed_left + 1 * crossed_down + 2 * crossed_right +
1064  3 * crossed_up;
1065  }
1066  else
1067  {
1068  incr = -1;
1069  insert = 3 * crossed_left + 2 * crossed_up + 1 * crossed_right +
1070  0 * crossed_down;
1071  }
1072  }
1073 
1074  if ( crossed_left + crossed_right == 2 && crossed_down + crossed_up == 0 )
1075  {
1076  if ( xclp[iclp - 1] == xmin )
1077  {
1078  if ( dir == 1 )
1079  {
1080  incr = 1;
1081  insert = 0;
1082  }
1083  else
1084  {
1085  incr = -1;
1086  insert = 3;
1087  }
1088  }
1089  else
1090  {
1091  if ( dir == 1 )
1092  {
1093  incr = 1;
1094  insert = 1;
1095  }
1096  else
1097  {
1098  incr = -1;
1099  insert = 2;
1100  }
1101  }
1102  }
1103 
1104  if ( crossed_left + crossed_right == 0 && crossed_down + crossed_up == 2 )
1105  {
1106  if ( yclp[iclp - 1] == ymin )
1107  {
1108  if ( dir == 1 )
1109  {
1110  incr = 1;
1111  insert = 1;
1112  }
1113  else
1114  {
1115  incr = -1;
1116  insert = 0;
1117  }
1118  }
1119  else
1120  {
1121  if ( dir == 1 )
1122  {
1123  incr = 1;
1124  insert = 3;
1125  }
1126  else
1127  {
1128  incr = -1;
1129  insert = 2;
1130  }
1131  }
1132  }
1133 
1134  for ( i = 0; i < 4; i++ )
1135  {
1136  xclp[iclp] = (short) xlim[insert];
1137  yclp[iclp] = (short) ylim[insert];
1138  iclp++;
1139  insert += incr;
1140  if ( insert > 3 )
1141  insert = 0;
1142  if ( insert < 0 )
1143  insert = 3;
1144  }
1145  }
1146 
1147  // Draw the sucker
1148  if ( iclp >= 3 )
1149  ( *draw )( xclp, yclp, iclp );
1150 
1151  if ( xclp != _xclp )
1152  {
1153  free( xclp );
1154  free( yclp );
1155  }
1156 }
1157 #endif // USE_FILL_INTERSECTION_POLYGON
1158 
1159 //--------------------------------------------------------------------------
1160 // int circulation()
1161 //
1162 // Returns the circulation direction for a given polyline: positive is
1163 // counterclockwise, negative is clockwise (right hand rule).
1164 //
1165 // Used to get the circulation of the fill polygon around the bounding box,
1166 // when the fill polygon is larger than the bounding box. Counts left
1167 // (positive) vs right (negative) hand turns using a cross product, instead of
1168 // performing all the expensive trig calculations needed to get this 100%
1169 // correct. For the fill cases encountered in plplot, this treatment should
1170 // give the correct answer most of the time, by far. When used with plshades,
1171 // the typical return value is 3 or -3, since 3 turns are necessary in order
1172 // to complete the fill region. Only for really oddly shaped fill regions
1173 // will it give the wrong answer.
1174 //
1175 // AM:
1176 // Changed the computation: use the outer product to compute the surface
1177 // area, the sign determines if the polygon is followed clockwise or
1178 // counterclockwise. This is more reliable. Floating-point numbers
1179 // are used to avoid overflow.
1180 //--------------------------------------------------------------------------
1181 
1182 int
1183 circulation( PLINT *x, PLINT *y, PLINT npts )
1184 {
1185  PLFLT xproduct;
1186  int direction = 0;
1187  PLFLT x1, y1, x2, y2, x3, y3;
1188  int i;
1189 
1190  xproduct = 0.0;
1191  x1 = x[0];
1192  y1 = y[0];
1193  for ( i = 1; i < npts - 2; i++ )
1194  {
1195  x2 = x[i + 1];
1196  y2 = y[i + 1];
1197  x3 = x[i + 2];
1198  y3 = y[i + 2];
1199  xproduct = xproduct + ( x2 - x1 ) * ( y3 - y2 ) - ( y2 - y1 ) * ( x3 - x2 );
1200  }
1201 
1202  if ( xproduct > 0.0 )
1203  direction = 1;
1204  if ( xproduct < 0.0 )
1205  direction = -1;
1206  return direction;
1207 }
1208 
1209 
1210 // PLFLT wrapper for !notpointinpolygon.
1211 int
1212 plP_pointinpolygon( PLINT n, const PLFLT *x, const PLFLT *y, PLFLT xp, PLFLT yp )
1213 {
1214  int i, return_value;
1215  PLINT *xint, *yint;
1216  PLFLT xmaximum = fabs( xp ), ymaximum = fabs( yp ), xscale, yscale;
1217  if ( ( xint = (PLINT *) malloc( (size_t) n * sizeof ( PLINT ) ) ) == NULL )
1218  {
1219  plexit( "PlP_pointinpolygon: Insufficient memory" );
1220  }
1221  if ( ( yint = (PLINT *) malloc( (size_t) n * sizeof ( PLINT ) ) ) == NULL )
1222  {
1223  plexit( "PlP_pointinpolygon: Insufficient memory" );
1224  }
1225  for ( i = 0; i < n; i++ )
1226  {
1227  xmaximum = MAX( xmaximum, fabs( x[i] ) );
1228  ymaximum = MAX( ymaximum, fabs( y[i] ) );
1229  }
1230  xscale = 1.e8 / xmaximum;
1231  yscale = 1.e8 / ymaximum;
1232  for ( i = 0; i < n; i++ )
1233  {
1234  xint[i] = (PLINT) ( xscale * x[i] );
1235  yint[i] = (PLINT) ( yscale * y[i] );
1236  }
1237  return_value = !notpointinpolygon( n, xint, yint,
1238  (PLINT) ( xscale * xp ), (PLINT) ( yscale * yp ) );
1239  free( xint );
1240  free( yint );
1241  return return_value;
1242 }
1243 //--------------------------------------------------------------------------
1244 // int notpointinpolygon()
1245 //
1246 // Returns 0, 1, or 2 depending on whether the test point is definitely
1247 // inside, near the border, or definitely outside the polygon.
1248 // Notes:
1249 // This "Ray casting algorithm" has been described in
1250 // http://en.wikipedia.org/wiki/Point_in_polygon.
1251 // Logic still needs to be inserted to take care of the "ray passes
1252 // through vertex" problem in a numerically robust way.
1253 //--------------------------------------------------------------------------
1254 
1255 // Temporary until get rid of old code altogether.
1256 #define NEW_NOTPOINTINPOLYGON_CODE
1257 static int
1258 notpointinpolygon( PLINT n, const PLINT *x, const PLINT *y, PLINT xp, PLINT yp )
1259 {
1260 #ifdef NEW_NOTPOINTINPOLYGON_CODE
1261  int i, im1, ifnotcrossed;
1262  int count_crossings = 0;
1263  PLINT xmin, xout, yout, xintersect, yintersect;
1264 
1265 
1266  // Determine a point outside the polygon
1267 
1268  xmin = x[0];
1269  xout = x[0];
1270  yout = y[0];
1271  for ( i = 1; i < n; i++ )
1272  {
1273  xout = MAX( xout, x[i] );
1274  xmin = MIN( xmin, x[i] );
1275  }
1276  // + 10 to make sure completely outside.
1277  xout = xout + ( xout - xmin ) + 10;
1278 
1279  // Determine whether the line between (xout, yout) and (xp, yp) intersects
1280  // one of the polygon segments.
1281 
1282  im1 = n - 1;
1283  for ( i = 0; i < n; i++ )
1284  {
1285  if ( !( x[im1] == x[i] && y[im1] == y[i] ) )
1286  {
1287  ifnotcrossed = notcrossed( &xintersect, &yintersect,
1288  x[im1], y[im1], x[i], y[i],
1289  xp, yp, xout, yout );
1290 
1291  if ( !ifnotcrossed )
1292  count_crossings++;
1293  else if ( ifnotcrossed & ( PL_NEAR_A1 | PL_NEAR_A2 | PL_NEAR_B1 | PL_NEAR_B2 ) )
1294  return 1;
1295  }
1296  im1 = i;
1297  }
1298 
1299  // return 0 if the test point is definitely inside
1300  // (count_crossings odd), return 1 if the test point is near (see
1301  // above logic), and return 2 if the test point is definitely
1302  // outside the border (count_crossings even).
1303  if ( ( count_crossings % 2 ) == 1 )
1304  return 0;
1305  else
1306  return 2;
1307 }
1308 #else // NEW_NOTPOINTINPOLYGON_CODE
1309  int i;
1310  int count_crossings;
1311  PLFLT x1, y1, x2, y2, xpp, ypp, xout, yout, xmax;
1312  PLFLT xvp, yvp, xvv, yvv, xv1, yv1, xv2, yv2;
1313  PLFLT inprod1, inprod2;
1314 
1315  xpp = (PLFLT) xp;
1316  ypp = (PLFLT) yp;
1317 
1318  count_crossings = 0;
1319 
1320 
1321  // Determine a point outside the polygon
1322 
1323  xmax = x[0];
1324  xout = x[0];
1325  yout = y[0];
1326  for ( i = 0; i < n; i++ )
1327  {
1328  if ( xout > x[i] )
1329  {
1330  xout = x[i];
1331  }
1332  if ( xmax < x[i] )
1333  {
1334  xmax = x[i];
1335  }
1336  }
1337  xout = xout - ( xmax - xout );
1338 
1339  // Determine for each side whether the line segment between
1340  // our two points crosses the vertex
1341 
1342  xpp = (PLFLT) xp;
1343  ypp = (PLFLT) yp;
1344 
1345  xvp = xpp - xout;
1346  yvp = ypp - yout;
1347 
1348  for ( i = 0; i < n; i++ )
1349  {
1350  x1 = (PLFLT) x[i];
1351  y1 = (PLFLT) y[i];
1352  if ( i < n - 1 )
1353  {
1354  x2 = (PLFLT) x[i + 1];
1355  y2 = (PLFLT) y[i + 1];
1356  }
1357  else
1358  {
1359  x2 = (PLFLT) x[0];
1360  y2 = (PLFLT) y[0];
1361  }
1362 
1363  // Skip zero-length segments
1364  if ( x1 == x2 && y1 == y2 )
1365  {
1366  continue;
1367  }
1368 
1369  // Line through the two fixed points:
1370  // Are x1 and x2 on either side?
1371  xv1 = x1 - xout;
1372  yv1 = y1 - yout;
1373  xv2 = x2 - xout;
1374  yv2 = y2 - yout;
1375  inprod1 = xv1 * yvp - yv1 * xvp; // Well, with the normal vector
1376  inprod2 = xv2 * yvp - yv2 * xvp;
1377  if ( inprod1 * inprod2 >= 0.0 )
1378  {
1379  // No crossing possible!
1380  continue;
1381  }
1382 
1383  // Line through the two vertices:
1384  // Are xout and xpp on either side?
1385  xvv = x2 - x1;
1386  yvv = y2 - y1;
1387  xv1 = xpp - x1;
1388  yv1 = ypp - y1;
1389  xv2 = xout - x1;
1390  yv2 = yout - y1;
1391  inprod1 = xv1 * yvv - yv1 * xvv;
1392  inprod2 = xv2 * yvv - yv2 * xvv;
1393  if ( inprod1 * inprod2 >= 0.0 )
1394  {
1395  // No crossing possible!
1396  continue;
1397  }
1398 
1399  // We do have a crossing
1400  count_crossings++;
1401  }
1402 
1403  // Return the result: an even number of crossings means the
1404  // point is outside the polygon
1405 
1406  return !( count_crossings % 2 );
1407 }
1408 #endif // NEW_NOTPOINTINPOLYGON_CODE
1409 
1410 #define MAX_RECURSION_DEPTH 10
1411 
1412 // Fill intersection of two simple polygons that do no self-intersect,
1413 // and which have no duplicate vertices or two consecutive edges that
1414 // are parallel. A further requirement is that both polygons have a
1415 // positive orientation (see
1416 // http://en.wikipedia.org/wiki/Curve_orientation). That is, as you
1417 // traverse the boundary in index order, the inside area of the
1418 // polygon is always on the left. Finally, the first vertex of
1419 // polygon 1 (starting with n1 -1) that is not near the border of
1420 // polygon 2 must be outside polygon 2. N.B. it is the calling
1421 // routine's responsibility to insure all those requirements are
1422 // satisfied.
1423 //
1424 // Two polygons that do not self intersect must have an even number of
1425 // edge crossings between them. (ignoring vertex intersections which
1426 // touch, but do not cross). fill_intersection_polygon eliminates
1427 // those intersection crossings by recursion (calling the same routine
1428 // twice again with the second polygon split at a boundary defined by
1429 // the first intersection point, all polygon 1 vertices between the
1430 // intersections, and the second intersection point). Once the
1431 // recursion has eliminated all crossing edges, fill or not using the
1432 // appropriate polygon depending on whether the first and second
1433 // polygons are identical or whether one of them is entirely inside
1434 // the other of them. If ifextrapolygon is true, the fill step will
1435 // consist of another recursive call to the routine with
1436 // ifextrapolygon false, and the second polygon set to an additional
1437 // polygon defined by the stream (not yet implemented).
1438 
1439 // arguments to intersection_polygon:
1440 // recursion_depth is just what it says.
1441 // ifextrapolygon used to decide whether to use extra polygon from the stream.
1442 //fill is the fill routine.
1443 //x1, *y1, n1 define the polygon 1 vertices.
1444 // i1start is the first polygon 1 index to look at (because all the previous
1445 // ones have been completely processed).
1446 //x2, *y2, *if2, n2 define the polygon 2 vertices plus a status indicator
1447 // for each vertex which is 1 for a previous crossing and 2 for a polygon
1448 // 1 vertex.
1449 // fill_status is 1 when polygons 1 and 2 _must_ include some joint
1450 // filled area and is -1 when polygons 1 and 2 _must_ include some
1451 // unfilled area. fill_status of +/- 1 is determined from the
1452 // orientations of polygon 1 and 2 from the next higher recursion
1453 // level and how those two are combined to form the polygon 2
1454 // split at this recursion level. fill_status = 0 occurs (at
1455 // recursion level 0) for polygons 1 and 2 that are independent of
1456 // each other.
1457 
1458 #ifdef USE_FILL_INTERSECTION_POLYGON
1459 void
1460 fill_intersection_polygon( PLINT recursion_depth, PLINT ifextrapolygon,
1461  PLINT fill_status,
1462  void ( *fill )( short *, short *, PLINT ),
1463  const PLINT *x1, const PLINT *y1,
1464  PLINT i1start, PLINT n1,
1465  const PLINT *x2, const PLINT *y2,
1466  const PLINT *if2, PLINT n2 )
1467 {
1468  PLINT i1, i1m1, i1start_new,
1469  i2, i2m1,
1470  kk, kkstart1, kkstart21, kkstart22,
1471  k, kstart, range1,
1472  range21, range22, ncrossed, ncrossed_change,
1473  nsplit1, nsplit2, nsplit2m1;
1474  PLINT xintersect[2], yintersect[2], i1intersect[2],
1475  i2intersect[2];
1476  PLINT *xsplit1, *ysplit1, *ifsplit1,
1477  *xsplit2, *ysplit2, *ifsplit2;
1478  PLINT ifill, nfill = 0,
1479  ifnotpolygon1inpolygon2, ifnotpolygon2inpolygon1;
1480  const PLINT *xfiller, *yfiller;
1481  short *xfill, *yfill;
1482 
1483  if ( recursion_depth > MAX_RECURSION_DEPTH )
1484  {
1485  plwarn( "fill_intersection_polygon: Recursion_depth too large. "
1486  "Probably an internal error figuring out intersections. " );
1487  return;
1488  }
1489 
1490  if ( n1 < 3 )
1491  {
1492  plwarn( "fill_intersection_polygon: Internal error; n1 < 3." );
1493  return;
1494  }
1495 
1496  if ( n2 < 3 )
1497  {
1498  plwarn( "fill_intersection_polygon: Internal error; n2 < 3." );
1499  return;
1500  }
1501 
1502  if ( i1start < 0 || i1start >= n1 )
1503  {
1504  plwarn( "fill_intersection_polygon: invalid i1start." );
1505  return;
1506  }
1507 
1508  // Check that there are no duplicate vertices.
1509  i1m1 = i1start - 1;
1510  if ( i1m1 < 0 )
1511  i1m1 = n1 - 1;
1512 
1513  for ( i1 = i1start; i1 < n1; i1++ )
1514  {
1515  if ( x1[i1] == x1[i1m1] && y1[i1] == y1[i1m1] )
1516  break;
1517  i1m1 = i1;
1518  }
1519 
1520  if ( i1 < n1 )
1521  {
1522  plwarn( "fill_intersection_polygon: Internal error; i1 < n1." );
1523  return;
1524  }
1525 
1526  i2m1 = n2 - 1;
1527  for ( i2 = 0; i2 < n2; i2++ )
1528  {
1529  if ( x2[i2] == x2[i2m1] && y2[i2] == y2[i2m1] )
1530  break;
1531  i2m1 = i2;
1532  }
1533 
1534  if ( i2 < n2 )
1535  {
1536  plwarn( "fill_intersection_polygon: Internal error; i2 < n2." );
1537  return;
1538  }
1539 
1540  //
1541  //
1542  // Follow polygon 1 (checking intersections with polygon 2 for each
1543  // segment of polygon 1) until you have accumulated two
1544  // intersections with polygon 2. Here is an ascii-art illustration
1545  // of the situation.
1546  //
1547  //
1548  // 2???2
1549  //
1550  // 2 2
1551  //
1552  // --- 1 1
1553  // 1 2 1 1 ...
1554  // X
1555  // 1
1556  // X
1557  // 2
1558  // 1 1
1559  // 1
1560  // 2
1561  // 2
1562  // 2???2
1563  //
1564  //
1565  // "1" marks polygon 1 vertices, "2" marks polygon 2 vertices, "X"
1566  // marks the intersections, "---" stands for part of polygon 1
1567  // that has been previously searched for all possible intersections
1568  // from index 0, and "..." means polygon 1 continues
1569  // with more potential intersections above and/or below this diagram
1570  // before it finally hooks back to connect with the index 0 vertex.
1571  // "2???2" stands for parts of polygon 2 that must connect with each other
1572  // (since the polygon 1 path between the two intersections is
1573  // known to be free of intersections.)
1574  //
1575  // Polygon 2 is split at the boundary defined by the two
1576  // intersections and all (in this case three) polygon 1 vertices
1577  // between the two intersections for the next recursion level. We
1578  // absolutely know for that boundary that no more intersections can
1579  // occur (both polygon 1 and polygon 2 are guaranteed not to
1580  // self-intersect) so we mark the status of those vertices with that
1581  // information so those polygon 2 split vertices will not be used to
1582  // search for further intersections at deeper recursion levels.
1583  // Note, we know nothing about whether the remaining "2???2" parts of the
1584  // split polygon 2 intersect with polygon 1 or not so those will
1585  // continued to be searched at deeper recursion levels. At the same
1586  // time, we absolutely know that the part of polygon 1 to the left of
1587  // rightmost x down to and including index 0 cannot yield more
1588  // intersections with any split of polygon 2 so we adjust the lower
1589  // limit of polygon 1 to be used for intersection searches at deeper
1590  // recursion levels. The result should be that at sufficiently deep
1591  // recursion depth we always end up with the case that there are no
1592  // intersections to be found between polygon 1 and some polygon 2
1593  // split, and in that case we move on to the end phase below.
1594  //
1595  ncrossed = 0;
1596  i1m1 = i1start - 1;
1597  if ( i1m1 < 0 )
1598  i1m1 += n1;
1599  for ( i1 = i1start; i1 < n1; i1++ )
1600  {
1601  ncrossed_change = number_crossings(
1602  &xintersect[ncrossed], &yintersect[ncrossed],
1603  &i2intersect[ncrossed], 2 - ncrossed,
1604  i1, n1, x1, y1, n2, x2, y2 );
1605  if ( ncrossed_change > 0 )
1606  {
1607  i1intersect[ncrossed] = i1;
1608  if ( ncrossed_change == 2 )
1609  {
1610  ;
1611  }
1612  i1intersect[1] = i1;
1613 
1614  ncrossed += ncrossed_change;
1615  if ( ncrossed == 2 )
1616  {
1617  // Have discovered the first two crossings for
1618  // polygon 1 at i1 = i1start or above.
1619 
1620  // New i1start is the same as the current i1 (just
1621  // in case there are more crossings to find between
1622  // i1m1 and i1.)
1623  i1start_new = i1;
1624 
1625  // Split polygon 2 at the boundary consisting of
1626  // first intersection, intervening (if any) range1
1627  // polygon 1 points and second intersection.
1628  // range1 must always be non-negative because i1
1629  // range only traversed once.
1630  range1 = i1intersect[1] - i1intersect[0];
1631  // Polygon 2 intersects could be anywhere (since
1632  // i2 range repeated until get an intersect).
1633  // Divide polygon 2 into two polygons with a
1634  // common boundary consisting of the first intersect,
1635  // range1 points from polygon 1 starting at index
1636  // kkstart1 of polygon 1, and the second intersect.
1637  kkstart1 = i1intersect[0];
1638 
1639  // Calculate polygon 2 index range in split 1 (the
1640  // split that proceeds beyond the second intersect with
1641  // ascending i2 values).
1642  range21 = i2intersect[0] - i2intersect[1];
1643  if ( range21 < 0 )
1644  range21 += n2;
1645  // i2 intersect values range between 0 and n2 - 1 so
1646  // the smallest untransformed range21 value is -n2 + 1,
1647  // and the largest untransformed range21 value is n2 - 1.
1648  // This means the smallest transformed range21 value is 0
1649  // (which occurs only ifi2intersect[0] = i2intersect[1],
1650  // see more commentary for that special case below) while
1651  // the largest transformed range21 value is n2 - 1.
1652 
1653  if ( range21 == 0 )
1654  {
1655  int ifxsort, ifascend;
1656  // For this case, the two crossings occur within the same
1657  // polygon 2 boundary segment and if those two crossings
1658  // are in ascending/descending order in i2, then split 1
1659  // (the split with the positive fill_status) must include
1660  // all/none of the points in polygon 2.
1661  i2 = i2intersect[1];
1662  i2m1 = i2 - 1;
1663  if ( i2m1 < 0 )
1664  i2m1 += n2;
1665 
1666  ifxsort = abs( x2[i2] - x2[i2m1] ) > abs( y2[i2] - y2[i2m1] );
1667  ifascend = ( ifxsort && x2[i2] > x2[i2m1] ) ||
1668  ( !ifxsort && y2[i2] > y2[i2m1] );
1669  if ( ( ifxsort && ifascend && xintersect[0] < xintersect[1] ) ||
1670  ( !ifxsort && ifascend && yintersect[0] < yintersect[1] ) ||
1671  ( ifxsort && !ifascend && xintersect[0] >= xintersect[1] ) ||
1672  ( !ifxsort && !ifascend && yintersect[0] >= yintersect[1] ) )
1673  {
1674  range21 = n2;
1675  }
1676  }
1677 
1678  kkstart21 = i2intersect[1];
1679  nsplit1 = 2 + range1 + range21;
1680 
1681  // Split 2 of polygon 2 consists of the
1682  // boundary + range22 (= n2 - range21) points
1683  // between kkstart22 (= i2intersect[1]-1) and i2intersect[0] in
1684  // descending order of polygon 2 indices.
1685  range22 = n2 - range21;
1686  // Starting i2 index of split 2.
1687  kkstart22 = i2intersect[1] - 1;
1688  if ( kkstart22 < 0 )
1689  kkstart22 += n2;
1690  nsplit2 = 2 + range1 + range22;
1691 
1692  if ( ( xsplit1 = (PLINT *) malloc( (size_t) nsplit1 * sizeof ( PLINT ) ) ) == NULL )
1693  {
1694  plexit( "fill_intersection_polygon: Insufficient memory" );
1695  }
1696  if ( ( ysplit1 = (PLINT *) malloc( (size_t) nsplit1 * sizeof ( PLINT ) ) ) == NULL )
1697  {
1698  plexit( "fill_intersection_polygon: Insufficient memory" );
1699  }
1700  if ( ( ifsplit1 = (PLINT *) malloc( (size_t) nsplit1 * sizeof ( PLINT ) ) ) == NULL )
1701  {
1702  plexit( "fill_intersection_polygon: Insufficient memory" );
1703  }
1704 
1705  if ( ( xsplit2 = (PLINT *) malloc( (size_t) nsplit2 * sizeof ( PLINT ) ) ) == NULL )
1706  {
1707  plexit( "fill_intersection_polygon: Insufficient memory" );
1708  }
1709  if ( ( ysplit2 = (PLINT *) malloc( (size_t) nsplit2 * sizeof ( PLINT ) ) ) == NULL )
1710  {
1711  plexit( "fill_intersection_polygon: Insufficient memory" );
1712  }
1713  if ( ( ifsplit2 = (PLINT *) malloc( (size_t) nsplit2 * sizeof ( PLINT ) ) ) == NULL )
1714  {
1715  plexit( "fill_intersection_polygon: Insufficient memory" );
1716  }
1717  // Common boundary between split1 and split2.
1718  // N.B. Although basic index arithmetic for
1719  // split 2 is done in negative orientation
1720  // order because the index is decrementing
1721  // relative to the index of split 2, actually
1722  // store results in reverse order to preserve
1723  // the positive orientation that by assumption
1724  // both polygon 1 and 2 have.
1725  k = 0;
1726  xsplit1[k] = xintersect[0];
1727  ysplit1[k] = yintersect[0];
1728  ifsplit1[k] = 1;
1729  nsplit2m1 = nsplit2 - 1;
1730  xsplit2[nsplit2m1 - k] = xintersect[0];
1731  ysplit2[nsplit2m1 - k] = yintersect[0];
1732  ifsplit2[nsplit2m1 - k] = 1;
1733  kstart = k + 1;
1734  kk = kkstart1;
1735  // No wrap checks on kk index below because
1736  // it must always be in valid range (since
1737  // polygon 1 traversed only once).
1738  for ( k = kstart; k < range1 + 1; k++ )
1739  {
1740  xsplit1[k] = x1[kk];
1741  ysplit1[k] = y1[kk];
1742  ifsplit1[k] = 2;
1743  xsplit2[nsplit2m1 - k] = x1[kk];
1744  ysplit2[nsplit2m1 - k] = y1[kk++];
1745  ifsplit2[nsplit2m1 - k] = 2;
1746  }
1747  xsplit1[k] = xintersect[1];
1748  ysplit1[k] = yintersect[1];
1749  ifsplit1[k] = 1;
1750  xsplit2[nsplit2m1 - k] = xintersect[1];
1751  ysplit2[nsplit2m1 - k] = yintersect[1];
1752  ifsplit2[nsplit2m1 - k] = 1;
1753 
1754  // Finish off collecting split1 using ascending kk
1755  // values.
1756  kstart = k + 1;
1757  kk = kkstart21;
1758  for ( k = kstart; k < nsplit1; k++ )
1759  {
1760  xsplit1[k] = x2[kk];
1761  ysplit1[k] = y2[kk];
1762  ifsplit1[k] = if2[kk++];
1763  if ( kk >= n2 )
1764  kk -= n2;
1765  }
1766 
1767  // N.B. the positive orientation of split1 is
1768  // preserved since the index order is the same
1769  // as that of polygon 2, and by assumption
1770  // that polygon and polygon 1 have identical
1771  // positive orientations.
1772  fill_intersection_polygon(
1773  recursion_depth + 1, ifextrapolygon, 1, fill,
1774  x1, y1, i1start_new, n1,
1775  xsplit1, ysplit1, ifsplit1, nsplit1 );
1776  free( xsplit1 );
1777  free( ysplit1 );
1778  free( ifsplit1 );
1779 
1780  // Finish off collecting split2 using descending kk
1781  // values.
1782  kk = kkstart22;
1783  for ( k = kstart; k < nsplit2; k++ )
1784  {
1785  xsplit2[nsplit2m1 - k] = x2[kk];
1786  ysplit2[nsplit2m1 - k] = y2[kk];
1787  ifsplit2[nsplit2m1 - k] = if2[kk--];
1788  if ( kk < 0 )
1789  kk += n2;
1790  }
1791 
1792  // N.B. the positive orientation of split2 is
1793  // preserved since the index order is the same
1794  // as that of polygon 2, and by assumption
1795  // that polygon and polygon 1 have identical
1796  // positive orientations.
1797  fill_intersection_polygon(
1798  recursion_depth + 1, ifextrapolygon, -1, fill,
1799  x1, y1, i1start_new, n1,
1800  xsplit2, ysplit2, ifsplit2, nsplit2 );
1801  free( xsplit2 );
1802  free( ysplit2 );
1803  free( ifsplit2 );
1804  return;
1805  }
1806  }
1807  i1m1 = i1;
1808  }
1809 
1810  if ( ncrossed != 0 )
1811  {
1812  plwarn( "fill_intersection_polygon: Internal error; ncrossed != 0." );
1813  return;
1814  }
1815 
1816  // This end phase is reached only if no crossings are found.
1817 
1818  // If a fill_status of +/- 1 is known, use that to fill or not since
1819  // +1 corresponds to all of polygon 2 inside polygon 1 and -1
1820  // corresponds to none of polygon 2 inside polygon 1.
1821  if ( fill_status == -1 )
1822  return;
1823  else if ( fill_status == 1 )
1824  {
1825  nfill = n2;
1826  xfiller = x2;
1827  yfiller = y2;
1828  }
1829  else if ( fill_status == 0 )
1830  //else if ( 1 )
1831  {
1832  if ( recursion_depth != 0 )
1833  {
1834  plwarn( "fill_intersection_polygon: Internal error; fill_status == 0 for recursion_depth > 0" );
1835  return;
1836  }
1837  // For this case (recursion level 0) the two polygons are
1838  // completely independent with no crossings between them or
1839  // edges constructed from one another.
1840  //
1841  // The intersection of polygon 2 and 1, must be either of them (in
1842  // which case fill with the inner one), or neither of them (in
1843  // which case don't fill at all).
1844 
1845  // Classify polygon 1 by looking for first vertex in polygon 1
1846  // that is definitely inside or outside polygon 2.
1847  for ( i1 = 0; i1 < n1; i1++ )
1848  {
1849  if ( ( ifnotpolygon1inpolygon2 =
1850  notpointinpolygon( n2, x2, y2, x1[i1], y1[i1] ) ) != 1 )
1851  break;
1852  }
1853 
1854  // Classify polygon 2 by looking for first vertex in polygon 2
1855  // that is definitely inside or outside polygon 1.
1856  ifnotpolygon2inpolygon1 = 1;
1857  for ( i2 = 0; i2 < n2; i2++ )
1858  {
1859  // Do not bother checking vertices already known to be on the
1860  // boundary with polygon 1.
1861  if ( !if2[i2] && ( ifnotpolygon2inpolygon1 =
1862  notpointinpolygon( n1, x1, y1, x2[i2], y2[i2] ) ) != 1 )
1863  break;
1864  }
1865 
1866  if ( ifnotpolygon2inpolygon1 == 0 && ifnotpolygon1inpolygon2 == 0 )
1867  plwarn( "fill_intersection_polygon: Internal error; no intersections found but each polygon definitely inside the other!" );
1868  else if ( ifnotpolygon2inpolygon1 == 2 && ifnotpolygon1inpolygon2 == 2 )
1869  // The polygons do not intersect each other so do not fill in this
1870  // case.
1871  return;
1872  else if ( ifnotpolygon2inpolygon1 == 0 )
1873  {
1874  // Polygon 2 definitely inside polygon 1.
1875  nfill = n2;
1876  xfiller = x2;
1877  yfiller = y2;
1878  }
1879  else if ( ifnotpolygon1inpolygon2 == 0 )
1880  {
1881  // Polygon 1 definitely inside polygon 2.
1882  nfill = n1;
1883  xfiller = x1;
1884  yfiller = y1;
1885  }
1886  else if ( ifnotpolygon2inpolygon1 == 1 && ifnotpolygon1inpolygon2 == 1 )
1887  {
1888  // Polygon 2 vertices near polygon 1 border and vice versa which
1889  // implies the polygons are identical.
1890  nfill = n2;
1891  xfiller = x2;
1892  yfiller = y2;
1893  }
1894  else
1895  {
1896  // Polygon 1 inscribed in polygon 2 or vice versa. This is normally
1897  // unlikely for two independent polygons so the implementation is
1898  // ToDo.
1899  plwarn( "fill_intersection_polygon: inscribed polygons are still ToDo" );
1900  }
1901  }
1902 
1903  if ( nfill > 0 )
1904  {
1905  if ( ( xfill = (short *) malloc( (size_t) nfill * sizeof ( short ) ) ) == NULL )
1906  {
1907  plexit( "fill_intersection_polygon: Insufficient memory" );
1908  }
1909  if ( ( yfill = (short *) malloc( (size_t) nfill * sizeof ( short ) ) ) == NULL )
1910  {
1911  plexit( "fill_intersection_polygon: Insufficient memory" );
1912  }
1913  for ( ifill = 0; ifill < nfill; ifill++ )
1914  {
1915  xfill[ifill] = (short) xfiller[ifill];
1916  yfill[ifill] = (short) yfiller[ifill];
1917  }
1918  ( *fill )( xfill, yfill, nfill );
1919  free( xfill );
1920  free( yfill );
1921  }
1922 
1923  return;
1924 }
1925 #endif
1926 
1927 // Returns a 0 status code
1928 // if the two line segments A, and B defined
1929 // by their end points (xA1, yA1, xA2, yA2, xB1, yB1, xB2, and yB2)
1930 // definitely (i.e., intersection point is not near the ends of either
1931 // of the line segments) cross each other. Otherwise, return a status
1932 // code specifying the type of non-crossing (i.e., completely
1933 // separate, near one of the ends, parallel).
1934 // Only if status = 0, return the actual
1935 // intersection via the argument list pointers to
1936 // xintersect and yintersect.
1937 
1938 int
1939 notcrossed( PLINT * xintersect, PLINT * yintersect,
1940  PLINT xA1, PLINT yA1, PLINT xA2, PLINT yA2,
1941  PLINT xB1, PLINT yB1, PLINT xB2, PLINT yB2 )
1942 {
1943  PLFLT factor, factor_NBCC, fxintersect, fyintersect;
1944  // These variables are PLFLT not for precision, but to
1945  // avoid integer overflows if they were typed as PLINT.
1946  PLFLT xA2A1, yA2A1, xB2B1, yB2B1;
1947  PLFLT xB1A1, yB1A1, xB2A1, yB2A1;
1948  PLINT status = 0;
1949  //
1950  // Two linear equations to be solved for x and y.
1951  // y = ((x - xA1)*yA2 + (xA2 - x)*yA1)/(xA2 - xA1)
1952  // y = ((x - xB1)*yB2 + (xB2 - x)*yB1)/(xB2 - xB1)
1953  //
1954  // Transform those two equations to coordinate system with origin
1955  // at (xA1, yA1).
1956  // y' = x'*yA2A1/xA2A1
1957  // y' = ((x' - xB1A1)*yB2A1 + (xB2A1 - x')*yB1A1)/xB2B1
1958  // ==>
1959  // x' = -(
1960  // (-xB1A1*yB2A1 + xB2A1*yB1A1)/xB2B1)/
1961  // (yB2B1/xB2B1 - yA2A1/xA2A1)
1962  // = (xB1A1*yB2A1 - xB2A1*yB1A1)*xA2A1/
1963  // (xA2A1*yB2B1 - yA2A1*xB2B1)
1964  //
1965  //
1966 
1967  xA2A1 = xA2 - xA1;
1968  yA2A1 = yA2 - yA1;
1969  xB2B1 = xB2 - xB1;
1970  yB2B1 = yB2 - yB1;
1971 
1972  factor = xA2A1 * yB2B1 - yA2A1 * xB2B1;
1973  factor_NBCC = PL_NBCC * ( fabs( xA2A1 ) + fabs( yB2B1 ) + fabs( yA2A1 ) + fabs( xB2B1 ) );
1974  if ( fabs( factor ) <= factor_NBCC )
1975  {
1976  if ( fabs( factor ) > 0. )
1977  status = status | PL_NEAR_PARALLEL;
1978  else
1979  status = status | PL_PARALLEL;
1980  }
1981  else
1982  {
1983  xB1A1 = xB1 - xA1;
1984  yB1A1 = yB1 - yA1;
1985  xB2A1 = xB2 - xA1;
1986  yB2A1 = yB2 - yA1;
1987 
1988  factor = ( xB1A1 * yB2A1 - yB1A1 * xB2A1 ) / factor;
1989  fxintersect = factor * xA2A1 + xA1;
1990  fyintersect = factor * yA2A1 + yA1;
1991  // The "redundant" x and y segment range checks (which include near the
1992  // end points) are needed for the vertical and the horizontal cases.
1993  if ( ( BETW_NBCC( fxintersect, xA1, xA2 ) && BETW_NBCC( fyintersect, yA1, yA2 ) ) &&
1994  ( BETW_NBCC( fxintersect, xB1, xB2 ) && BETW_NBCC( fyintersect, yB1, yB2 ) ) )
1995  {
1996  // The intersect is close (within +/- PL_NBCC) to an end point or
1997  // corresponds to a definite crossing of the two line segments.
1998  // Find out which.
1999  if ( fabs( fxintersect - xA1 ) <= PL_NBCC && fabs( fyintersect - yA1 ) <= PL_NBCC )
2000  status = status | PL_NEAR_A1;
2001  else if ( fabs( fxintersect - xA2 ) <= PL_NBCC && fabs( fyintersect - yA2 ) <= PL_NBCC )
2002  status = status | PL_NEAR_A2;
2003  else if ( fabs( fxintersect - xB1 ) <= PL_NBCC && fabs( fyintersect - yB1 ) <= PL_NBCC )
2004  status = status | PL_NEAR_B1;
2005  else if ( fabs( fxintersect - xB2 ) <= PL_NBCC && fabs( fyintersect - yB2 ) <= PL_NBCC )
2006  status = status | PL_NEAR_B2;
2007  // N.B. if none of the above conditions hold then status remains at
2008  // zero to signal we have a definite crossing.
2009  }
2010  else
2011  status = status | PL_NOT_CROSSED;
2012  }
2013  if ( !status )
2014  {
2015  *xintersect = (PLINT) fxintersect;
2016  *yintersect = (PLINT) fyintersect;
2017  }
2018 
2019  return status;
2020 }
2021 
2022 #ifdef USE_FILL_INTERSECTION_POLYGON
2023 // Decide if polygon has a positive orientation or not.
2024 // Note the simple algorithm given in
2025 // http://en.wikipedia.org/wiki/Curve_orientation is incorrect for
2026 // non-convex polygons. Instead, for the general nonintersecting case
2027 // use the polygon area method given by
2028 // http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/ where
2029 // you project each edge of the polygon down to the X axis and take the
2030 // area of the enclosed trapezoid. The trapezoid areas outside the
2031 // polygon are completely cancelled if you sum over all edges. Furthermore,
2032 // the sum of the trapezoid areas have terms which are zero when calculated
2033 // with the telescoping rule so the final formula is quite simple.
2034 int
2035 positive_orientation( PLINT n, const PLINT *x, const PLINT *y )
2036 {
2037  PLINT i, im1;
2038  // Use PLFLT for all calculations to avoid integer overflows. Also,
2039  // the normal PLFLT has 64-bits which means you get exact integer
2040  // arithmetic well beyond the normal integer overflow limits.
2041  PLFLT twice_area = 0.;
2042  if ( n < 3 )
2043  {
2044  plwarn( "positive_orientation: internal logic error, n < 3" );
2045  return 0;
2046  }
2047  im1 = n - 1;
2048  for ( i = 0; i < n; i++ )
2049  {
2050  twice_area += (PLFLT) x[im1] * (PLFLT) y[i] - (PLFLT) x[i] * (PLFLT) y[im1];
2051  im1 = i;
2052  }
2053  if ( twice_area == 0. )
2054  {
2055  plwarn( "positive_orientation: internal logic error, twice_area == 0." );
2056  return 0;
2057  }
2058  else
2059  return twice_area > 0.;
2060 }
2061 
2062 // For the line with endpoints which are the (i1-1)th, and i1th
2063 // vertices of polygon 1 with polygon 2 find all definite crossings
2064 // with polygon 1. (The full polygon 1 information is needed only to
2065 // help sort out (NOT IMPLEMENTED YET) ambiguous crossings near a
2066 // vertex of polygon 1.) Sort those definite crossings in ascending
2067 // order along the line between the (i1-1)th and i1th vertices of
2068 // polygon 1, and return the first ncross (= 1 or = 2) crossings in the
2069 // xcross, ycross, and i2cross arrays. Return a zero or positive
2070 // status code of the actual number of crossings retained up to the
2071 // maximum allowed value of ncross. If some error occurred, return a
2072 // negative status code.
2073 
2074 int
2075 number_crossings( PLINT *xcross, PLINT *ycross, PLINT *i2cross, PLINT ncross,
2076  PLINT i1, PLINT n1, const PLINT *x1, const PLINT *y1,
2077  PLINT n2, const PLINT *x2, const PLINT *y2 )
2078 {
2079  int i1m1, i2, i2m1, ifnotcrossed;
2080  int ifxsort, ifascend, count_crossings = 0, status = 0;
2081  PLINT xintersect, yintersect;
2082 
2083  i1m1 = i1 - 1;
2084  if ( i1m1 < 0 )
2085  i1m1 += n1;
2086  if ( !( ncross == 1 || ncross == 2 ) ||
2087  ( x1[i1m1] == x1[i1] && y1[i1m1] == y1[i1] ) || n1 < 2 || n2 < 2 )
2088  {
2089  plwarn( "findcrossings: invalid call" );
2090  return -1;
2091  }
2092 
2093  ifxsort = abs( x1[i1] - x1[i1m1] ) > abs( y1[i1] - y1[i1m1] );
2094  ifascend = ( ifxsort && x1[i1] > x1[i1m1] ) ||
2095  ( !ifxsort && y1[i1] > y1[i1m1] );
2096 
2097  // Determine all crossings between the line between the (i1-1)th
2098  // and i1th vertices of polygon 1 and all edges of polygon 2.
2099  // Retain the lowest and (if ncross ==2) next lowest crossings in
2100  // order of x (or y if ifxsort is false) along the line from i1-1
2101  // to i1.
2102 
2103  i1m1 = i1 - 1;
2104  if ( i1m1 < 0 )
2105  i1m1 += n1;
2106  i2m1 = n2 - 1;
2107  for ( i2 = 0; i2 < n2; i2++ )
2108  {
2109  if ( !( x2[i2m1] == x2[i2] && y2[i2m1] == y2[i2] ) )
2110  {
2111  ifnotcrossed = notcrossed( &xintersect, &yintersect,
2112  x1[i1m1], y1[i1m1], x1[i1], y1[i1],
2113  x2[i2m1], y2[i2m1], x2[i2], y2[i2] );
2114 
2115  if ( !ifnotcrossed )
2116  {
2117  count_crossings++;
2118  if ( count_crossings == 1 )
2119  {
2120  xcross[0] = xintersect;
2121  ycross[0] = yintersect;
2122  i2cross[0] = i2;
2123  status = 1;
2124  }
2125  else
2126  {
2127  if ( ( ifxsort && ifascend && xintersect < xcross[0] ) ||
2128  ( !ifxsort && ifascend && yintersect < ycross[0] ) ||
2129  ( ifxsort && !ifascend && xintersect >= xcross[0] ) ||
2130  ( !ifxsort && !ifascend && yintersect >= ycross[0] ) )
2131  {
2132  if ( ncross == 2 )
2133  {
2134  xcross[1] = xcross[0];
2135  ycross[1] = ycross[0];
2136  i2cross[1] = i2cross[0];
2137  status = 2;
2138  }
2139  xcross[0] = xintersect;
2140  ycross[0] = yintersect;
2141  i2cross[0] = i2;
2142  }
2143  else if ( ncross == 2 && ( count_crossings == 2 || (
2144  ( ifxsort && ifascend && xintersect < xcross[1] ) ||
2145  ( !ifxsort && ifascend && yintersect < ycross[1] ) ||
2146  ( ifxsort && !ifascend && xintersect >= xcross[1] ) ||
2147  ( !ifxsort && !ifascend && yintersect >= ycross[1] ) ) ) )
2148  {
2149  xcross[1] = xintersect;
2150  ycross[1] = yintersect;
2151  i2cross[1] = i2;
2152  status = 2;
2153  }
2154  }
2155  }
2156  }
2157  i2m1 = i2;
2158  }
2159  return status;
2160 }
2161 #endif