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