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