/************************************************************************/
/*                */
/*  Copyright Pittsburgh Supercomputing Center 1987, 1988, 1990  */
/*  All Rights Reserved            */
/*  Author Joel Welling            */
/*                    */
/*      cgm_clip.c version 1.0                                          */
/*                                                                      */
/************************************************************************/

# include <stdlib.h>
# include <stdio.h>
# include "cgm_clip.h"

/* Flag to indicate whether or not clipping is to be performed */
static int clip_flag= 1;

/* Coordinate buffers and run length buffer, enlarged as needed by getcmem */
static int *ixbuf, *iybuf, *lengths;
static int ixbuf_sz=0, iybuf_sz=0, lengths_sz= 0;

/* Debugging control */
static int debug= 0;

/*
This macro returns true if the point (qx,qy) lies in the left half-plane
of the plane divided by the vector from (px,py) to (vx,vy).
*/
#define leftside( px, py, vx, vy, qx, qy ) \
  ( 0.0 <= ( (qx-px)*(-vy+py) + (qy-py)*(vx-px) ) )

/*
This macro sets the parameters npx and npy to the x and y coordinates
of the intersection point between a line from (px,py) to (qx,qy) and
a line at constant x=xval.
*/
#define y_intersect( px, py, qx, qy, xval, npx, npy ) \
  { npx= xval; npy= py + ((xval-px)*(qy-py)/(qx-px)); }

/*
This macro sets the parameters npx and npy to the x and y coordinates
of the intersection point between a line from (px,py) to (qx,qy) and
a line at constant y=yval.
*/
#define x_intersect( px, py, qx, qy, yval, npx, npy ) \
  { npy= yval; npx= px + ((yval-py)*(qx-px)/(qy-py)); }

void clip_on()
/* This function turns clipping on */
{
  if (debug) fprintf(stderr,"clip_on: clipping turned on.\n");
  clip_flag= 1;
}

void clip_off()
/* This function turns clipping off */
{
  if (debug) fprintf(stderr,"clip_off: clipping turned off.\n");
  clip_flag= 0;
}

void clip_toggle_debug()
/* This function turns clipping debugging on if it's off, or vice versa. */
{
  if (debug) {
    fprintf(stderr,"clip_toggle_debug: debugging off\n");
    debug= 0;
  }
  else {
    fprintf(stderr,"clip_toggle_debug: debugging on\n");
    debug= 1;
  }
}

static int getcmem(nx,ny,nlengths)
int nx,ny,nlengths;
/*
  This function checks that sufficient memory is available for
  coordinate storage, allocating more as needed.  Returns 1 if
  successful.
*/  
{
  if (nx>ixbuf_sz)
    {
    if (debug) fprintf(stderr,
       " getcmem: increasing x coord mem to %d ints.\n",nx);
    if (ixbuf_sz>0) { free((char *)ixbuf); ixbuf_sz= 0; };
    if ( (ixbuf= (int *)malloc(nx*sizeof(int))) == 0 )
      { 
      fprintf(stderr,"Unable to allocate x memory.\n");
       return(2);
      };
    ixbuf_sz= nx;
    };
  if (ny>iybuf_sz)
    {
    if (debug) fprintf(stderr,
       " getcmem: increasing y coord mem to %d ints.\n",ny);
    if (iybuf_sz>0) { free((char *)iybuf); iybuf_sz= 0; };
    if ( (iybuf= (int *)malloc(ny*sizeof(int))) == 0 )
      { 
      fprintf(stderr,"Unable to allocate y memory.\n");
      return(2);
      };
    iybuf_sz= ny;
    };
  if (nlengths>lengths_sz)
    {
    if (debug) fprintf(stderr,
    " getcmem: increasing run length buffer mem to %d ints.\n",nlengths);
    if (lengths_sz>0) { free((char *)lengths); lengths_sz= 0; };
    if ( (lengths= (int *)malloc(nlengths*sizeof(int))) == 0 )
      { 
      fprintf(stderr,"Unable to allocate length memory.\n");
      return(2);
      };
    lengths_sz= nlengths;
    };
  return(1);
}

static int no_clip( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
         xfun, yfun )
int npts, (*xfun)(), (*yfun)();
float *xarray, *yarray, xmin, xmax, ymin, ymax;
/*
This routine simply translates the coordinates without clipping.
*/
{
  int ipt, *ixcpy, *iycpy;

  if (debug) fprintf(stderr, "no_clip:\n");

  ixcpy= ixbuf; 
  iycpy= iybuf;
  for ( ipt= 0; ipt<npts; ipt++ ) {
    *ixcpy++= (*xfun)( *xarray++ );
    *iycpy++= (*yfun)( *yarray++ );
  };

  return(npts);
}

int clip_markers( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
      xfun, yfun, pixbuf, piybuf )
int npts, **pixbuf, **piybuf, (*xfun)(), (*yfun)();
float *xarray, *yarray, xmin, xmax, ymin, ymax;
/*
This routine converts a coordinate list to device (integer) coordinates,
clipping as for a polymarker.  The value returned is the total number
of markers to be plotted.
*/
{
  int ipt, *ixcpy, *iycpy, count;
  
  if (debug) fprintf(stderr, "clip_markers:\n");

  if ( getcmem(npts,npts,0) != 1 ) {
    fprintf(stderr," Error allocating memory in WRTPMK; call ignored.\n");
    return( 0 );
  };

  ixcpy= ixbuf; 
  iycpy= iybuf;
  if (clip_flag) {
    for ( ipt= 0; ipt<npts; ipt++ ) {
      if ( (*xarray >= xmin) && (*xarray <= xmax)
    && (*yarray >= ymin) && (*yarray <= ymax) ) {
  *ixcpy++= (*xfun)( *xarray );
  *iycpy++= (*yfun)( *yarray );
      }
      xarray++;
      yarray++;
    }
    count= ixcpy - ixbuf;
  }
  else 
    count= no_clip( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
         xfun, yfun );

  *pixbuf= ixbuf;
  *piybuf= iybuf;
  return( count );
}

static int clip_one_line( px, py, qx, qy, xmin, xmax, ymin, ymax,
       npx, npy, nqx, nqy, used_p, used_q, idcode )
float px, py, qx, qy, xmin, xmax, ymin, ymax;
float *npx, *npy, *nqx, *nqy;
int *used_p, *used_q, *idcode;
/*
This routine performs clipping on one line segment.  The line running
from (px, py) to (qx, qy) is clipped to the rectangle defined by
the x range xmin to xmax and the y range ymin to ymax.  New p and q
coordinates are returned via the pointers npx, npy, nqx, and nqy.  If
the new p coordinates are the same as the originals (i.e. the point
p lies within the clipping region) *used_p is set to non-zero;
likewise *used_q is non-zero if q was not clipped.  The function returns
1 if any part of the line appears within the clipping boundary, 0
otherwise.  If the function returns zero, used_p and used_q are set
but its other output parameters are not set.

The algorithm used is a cross between the Cohen-Sutherland and
Nicholl-Lee-Nicholl algorithms (see Computer Graphics, Principles
and Practice by Foley, van Dam, Feiner, and Hughes).  Points p and
q are classed according to which of the 9 regions relative to the
clip boundaries they lie in.  This makes 9x9=81 possible cases.
For each case, we check which edge to clip against by seeing which
side of a line from p to some appropriate vertex the point q lies
on, then clip.
                                    1   2   3
The regions are numbered like:      4   5   6
                                    7   8   9

where region 5 is within the clipping boundaries.
*/
{
  float tpx, tpy, tqx, tqy;

  /* Establish which case this is */
  *idcode= 0;

  if ( qy > ymax ) *idcode += 8;
  else if ( qy < ymin ) *idcode += 4;

  if ( qx > xmax ) *idcode += 2;
  else if (qx < xmin ) *idcode += 1;

  if ( py > ymax ) *idcode += 128;
  else if ( py < ymin ) *idcode += 64;

  if ( px > xmax ) *idcode += 32;
  else if ( px < xmin ) *idcode += 16;

  if (*idcode & 240) /* p not in region 5 */
    *used_p= 0;
  else *used_p= 1;

  if (*idcode & 15) /* q not in region 5 */
    *used_q= 0;
  else *used_q=1;

  /* The following macros are intended to be used in this routine only */
#define leftTL leftside( px, py, xmin, ymax, qx, qy )
#define leftTR leftside( px, py, xmax, ymax, qx, qy )
#define leftBL leftside( px, py, xmin, ymin, qx, qy )
#define leftBR leftside( px, py, xmax, ymin, qx, qy )
#define clipT( x, y ) x_intersect( px, py, qx, qy, ymax, x, y )
#define clipB( x, y ) x_intersect( px, py, qx, qy, ymin, x, y )
#define clipL( x, y ) y_intersect( px, py, qx, qy, xmin, x, y )
#define clipR( x, y ) y_intersect( px, py, qx, qy, xmax, x, y )

  switch (*idcode) {

  /* p in 1 */
  case 0x99: /* q in 1 */ return(0);
  case 0x98: /* q in 2 */ return(0);
  case 0x9a: /* q in 3 */ return(0);
  case 0x91: /* q in 4 */ return(0);
  case 0x90: /* q in 5 */ 
    *nqx= qx; *nqy= qy;
    if ( leftTL ) clipT( *npx, *npy )
    else clipL( *npx, *npy );
    return(1);
  case 0x92: /* q in 6 */ 
    if ( leftTR ) return(0);
    clipR( *nqx, *nqy );
    if ( leftTL ) clipT( *npx, *npy )
    else clipL( *npx, *npy );
    return(1);
  case 0x95: /* q in 7 */ return(0);
  case 0x94: /* q in 8 */ 
    if ( !leftBL ) return(0);
    clipB( *nqx, *nqy );
    if ( leftTL ) clipT( *npx, *npy )
    else clipL( *npx, *npy );
    return(1);
  case 0x96: /* q in 9 */ 
    if ( leftTR ) return(0);
    if ( !leftBL ) return(0);
    if ( leftTL ) {
      clipT( *npx, *npy );
      if ( leftBR ) clipR( *nqx, *nqy )
      else clipB( *nqx, *nqy );
    }
    else {
      clipL( *npx, *npy );
      if ( leftBR ) clipR( *nqx, *nqy )
      else clipB( *nqx, *nqy );
    }
    return(1);

  /* p in 2 */
  case 0x89: /* q in 1 */ return(0);
  case 0x88: /* q in 2 */ return(0);
  case 0x8a: /* q in 3 */ return(0);
  case 0x81: /* q in 4 */ 
    if ( !leftTL ) return(0);
    clipT( *npx, *npy );
    clipL( *nqx, *nqy );
    return(1);
  case 0x80: /* q in 5 */ 
    clipT( *npx, *npy );
    *nqx= qx; *nqy= qy;
    return(1);
  case 0x82: /* q in 6 */ 
    if ( leftTR ) return(0);
    clipT( *npx, *npy );
    clipR( *nqx, *nqy );
    return(1);
  case 0x85: /* q in 7 */ 
    if ( !leftTL ) return(0);
    clipT( *npx, *npy );
    if ( leftBL ) clipB( *nqx, *nqy )
    else clipL( *nqx, *nqy );
    return(1);
  case 0x84: /* q in 8 */ 
    clipT( *npx, *npy );
    clipB( *nqx, *nqy );
    return(1);
  case 0x86: /* q in 9 */
    if ( leftTR ) return(0);
    clipT( *npx, *npy );
    if ( leftBR ) clipR( *nqx, *nqy )
    else clipB( *nqx, *nqy );
    return(1);

  /* p in 3 */
  case 0xa9: /* q in 1 */ return(0);
  case 0xa8: /* q in 2 */ return(0);
  case 0xaa: /* q in 3 */ return(0);
  case 0xa1: /* q in 4 */ 
    if ( !leftTL ) return(0);
    clipL( *nqx, *nqy );
    if ( leftTR ) clipR( *npx, *npy )
    else clipT( *npx, *npy );
    return(1);
  case 0xa0: /* q in 5 */ 
    *nqx= qx; *nqy= qy;
    if ( leftTR ) clipR( *npx, *npy )
    else clipT( *npx, *npy );
    return(1);
  case 0xa2: /* q in 6 */ return(0);
  case 0xa5: /* q in 7 */ 
    if ( !leftTL ) return(0);
    if ( leftBR ) return(0);
    if ( leftTR ) {
      clipR( *npx, *npy );
      if ( leftBL ) clipB( *nqx, *nqy )
      else clipL( *nqx, *nqy );
    }
    else {
      clipT( *npx, *npy );
      if ( leftBL ) clipB( *nqx, *nqy )
      else clipL( *nqx, *nqy );
    }
    return(1);
  case 0xa4: /* q in 8 */ 
    if ( leftBR ) return(0);
    clipB( *nqx, *nqy );
    if ( leftTR ) clipR( *npx, *npy )
    else clipT( *npx, *npy );
    return(1);
  case 0xa6: /* q in 9 */ return(0);

  /* p in 4 */
  case 0x19: /* q in 1 */ return(0);
  case 0x18: /* q in 2 */ 
    if ( leftTL ) return(0);
    clipL( *npx, *npy );
    clipT( *nqx, *nqy );
    return(1);
  case 0x1a: /* q in 3 */ 
    if ( leftTL ) return(0);
    clipL( *npx, *npy );
    if ( leftTR ) clipT( *nqx, *nqy )
    else clipR( *nqx, *nqy );
    return(1);
  case 0x11: /* q in 4 */ return(0);
  case 0x10: /* q in 5 */ 
    clipL( *npx, *npy );
    *nqx= qx; *nqy= qy;
    return(1);
  case 0x12: /* q in 6 */ 
    clipL( *npx, *npy );
    clipR( *nqx, *nqy );
    return(1);
  case 0x15: /* q in 7 */ return(0);
  case 0x14: /* q in 8 */ 
    if ( !leftBL ) return(0);
    clipL( *npx, *npy );
    clipB( *nqx, *nqy );
    return(1);
  case 0x16: /* q in 9 */ 
    if ( !leftBL ) return(0);
    clipL( *npx, *npy );
    if ( leftBR ) clipR( *nqx, *nqy )
    else clipB( *nqx, *nqy );
    return(1);

  /* p in 5 */
  case 0x09: /* q in 1 */ 
    *npx= px; *npy= py;
    if ( leftTL ) clipL( *nqx, *nqy )
    else clipT( *nqx, *nqy )
    return(1);
  case 0x08: /* q in 2 */ 
    *npx= px; *npy= py;
    clipT( *nqx, *nqy );
    return(1);
  case 0x0a: /* q in 3 */ 
    *npx= px; *npy= py;
    if ( leftTR ) clipT( *nqx, *nqy )
    else clipR( *nqx, *nqy )
    return(1);
  case 0x01: /* q in 4 */ 
    *npx= px; *npy= py;
    clipL( *nqx, *nqy );
    return(1);
  case 0x00: /* q in 5 - case of trivial acceptance */ 
    *npx= px; *npy= py;
    *nqx= qx; *nqy= qy;
    return(1);
  case 0x02: /* q in 6 */ 
    *npx= px; *npy= py;
    clipR( *nqx, *nqy );
    return(1);
  case 0x05: /* q in 7 */ 
    *npx= px; *npy= py;
    if ( leftBL ) clipB( *nqx, *nqy )
    else clipL( *nqx, *nqy )
    return(1);
  case 0x04: /* q in 8 */ 
    *npx= px; *npy= py;
    clipB( *nqx, *nqy );
    return(1);
  case 0x06: /* q in 9 */ 
    *npx= px; *npy= py;
    if ( leftBR ) clipR( *nqx, *nqy )
    else clipB( *nqx, *nqy )
    return(1);

  /* p in 6 */
  case 0x29: /* q in 1 */ 
    if ( !leftTR ) return(0);
    clipR( *npx, *npy );
    if ( leftTL ) clipL( *nqx, *nqy )
    else clipT( *nqx, *nqy );
    return(1);
  case 0x28: /* q in 2 */ 
    if ( !leftTR ) return(0);
    clipR( *npx, *npy );
    clipT( *nqx, *nqy );
    return(1);
  case 0x2a: /* q in 3 */ return(0);
  case 0x21: /* q in 4 */ 
    clipR( *npx, *npy );
    clipL( *nqx, *nqy );
    return(1);
  case 0x20: /* q in 5 */ 
    clipR( *npx, *npy );
    *nqx= qx; *nqy= qy;
    return(1);
  case 0x22: /* q in 6 */ return(0);
  case 0x25: /* q in 7 */ 
    if ( leftBR ) return(0);
    clipR( *npx, *npy );
    if ( leftBL ) clipB( *nqx, *nqy )
    else clipL( *nqx, *nqy );
    return(1);
  case 0x24: /* q in 8 */ 
    if ( leftBR ) return(0);
    clipR( *npx, *npy );
    clipB( *nqx, *nqy );
    return(1);
  case 0x26: /* q in 9 */ return(0);

  /* p in 7 */
  case 0x59: /* q in 1 */ return(0);
  case 0x58: /* q in 2 */ 
    if ( leftTL ) return(0);
    clipT( *nqx, *nqy );
    if ( leftBL ) clipL( *npx, *npy )
    else clipB( *npx, *npy );
    return(1);
  case 0x5a: /* q in 3 */ 
    if ( leftTL ) return(0);
    if ( !leftBR ) return(0);
    if ( leftBL ) {
      clipL( *npx, *npy );
      if ( leftTR ) clipT( *nqx, *nqy )
      else clipR( *nqx, *nqy );
    }
    else {
      clipB( *npx, *npy );
      if ( leftTR ) clipT( *nqx, *nqy )
      else clipR( *nqx, *nqy );
    }
    return(1);
  case 0x51: /* q in 4 */ return(0);
  case 0x50: /* q in 5 */ 
    *nqx= qx; *nqy= qy;
    if ( leftBL ) clipL( *npx, *npy )
    else clipB( *npx, *npy );
    return(1);
  case 0x52: /* q in 6 */ 
    if ( !leftBR ) return(0);
    clipR( *nqx, *nqy );
    if ( leftBL ) clipL( *npx, *npy )
    else clipB( *npx, *npy );
    return(1);
  case 0x55: /* q in 7 */ return(0);
  case 0x54: /* q in 8 */ return(0);
  case 0x56: /* q in 9 */ return(0);

  /* p in 8 */
  case 0x49: /* q in 1 */ 
    if ( leftBL ) return(0);
    clipB( *npx, *npy );
    if ( leftTL ) clipL( *nqx, *nqy )
    else clipT( *nqx, *nqy );
    return(1);
  case 0x48: /* q in 2 */ 
    clipB( *npx, *npy );
    clipT( *nqx, *nqy );
    return(1);
  case 0x4a: /* q in 3 */ 
    if ( !leftBR ) return(0);
    clipB( *npx, *npy );
    if ( leftTR ) clipT( *nqx, *nqy )
    else clipR( *nqx, *nqy );
    return(1);
  case 0x41: /* q in 4 */ 
    if ( leftBL ) return(0);
    clipB( *npx, *npy );
    clipL( *nqx, *nqy );
    return(1);
  case 0x40: /* q in 5 */ 
    clipB( *npx, *npy );
    *nqx= qx; *nqy= qy;
    return(1);
  case 0x42: /* q in 6 */ 
    if ( !leftBR ) return(0);
    clipB( *npx, *npy );
    clipR( *nqx, *nqy );
    return(1);
  case 0x45: /* q in 7 */ return(0);
  case 0x44: /* q in 8 */ return(0);
  case 0x46: /* q in 9 */ return(0);

  /* p in 9 */
  case 0x69: /* q in 1 */ 
    if ( leftBL ) return(0);
    if ( !leftTR ) return(0);
    if ( leftBR ) {
      clipB( *npx, *npy );
      if ( leftTL ) clipL( *nqx, *nqy )
      else clipT( *nqx, *nqy );
    }
    else {
      clipR( *npx, *npy );
      if ( leftTL ) clipL( *nqx, *nqy )
      else clipT( *nqx, *nqy );
    }
    return(1);
  case 0x68: /* q in 2 */ 
    if ( !leftTR ) return(0);
    clipT( *nqx, *nqy );
    if ( leftBR ) clipB( *npx, *npy )
    else clipR( *npx, *npy );
    return(1);
  case 0x6a: /* q in 3 */ return(0);
  case 0x61: /* q in 4 */ 
    if ( leftBL ) return(0);
    clipL( *nqx, *nqy );
    if ( leftBR ) clipB( *npx, *npy )
    else clipR( *npx, *npy );
    return(1);
  case 0x60: /* q in 5 */ 
    *nqx= qx; *nqy= qy;
    if ( leftBR ) clipB( *npx, *npy )
    else clipR( *npx, *npy );
    return(1);
  case 0x62: /* q in 6 */ return(0);
  case 0x65: /* q in 7 */ return(0);
  case 0x64: /* q in 8 */ return(0);
  case 0x66: /* q in 9 */ return(0);

  default:
    fprintf(stderr,"invalid code 0x%x in clip_one_line!\n",*idcode);
    exit(2);
  }

  /* Undo the macro definitions */
#undef leftTL
#undef leftTR
#undef leftBL
#undef leftBR
#undef clipT
#undef clipB
#undef clipL
#undef clipR
  /*NOTREACHED*/
}

int clip_lines( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
      xfun, yfun, pixbuf, piybuf, plengths )
int npts, **pixbuf, **piybuf, **plengths, (*xfun)(), (*yfun)();
float *xarray, *yarray, xmin, xmax, ymin, ymax;
/*
This routine converts a coordinate list to device (integer) coordinates,
clipping as for a polyline.  The value returned is the total number
of polylines to be plotted;  the number of vertices in each polyline
is returned in the array pointed to by *plengths.
*/
{
  int ipt, *ixcpy, *iycpy, *lengthscpy, count;
  int used_p, used_q, idcode, run_length;
  float px, py, qx, qy, npx, npy, nqx, nqy;
  
  if (debug) fprintf(stderr, "clip_lines:\n");

  if (npts < 2) {
    fprintf(stderr," clip_lines: got only %d vertices, need at least 2!\n",
      npts);
    return( 0 );
  }

  if ( getcmem( 2*npts, 2*npts, npts ) != 1 ) {
    fprintf(stderr," Error allocating memory in WRTPMK; call ignored.\n");
    return( -1 );
  };

  ixcpy= ixbuf; 
  iycpy= iybuf;
  lengthscpy= lengths;
  if (clip_flag) {

    /* initialization */
    px= *xarray++;
    py= *yarray++;
    qx= *xarray++;
    qy= *yarray++;
    run_length= 1;
    count= 0;

    /* first two points */
    if ( clip_one_line( px, py, qx, qy, xmin, xmax, ymin, ymax,
      &npx, &npy, &nqx, &nqy, &used_p, &used_q, &idcode ) ) {
      run_length++;
      *ixcpy++= (*xfun)(npx); *iycpy++= (*yfun)(npy);
      *ixcpy++= (*xfun)(nqx); *iycpy++= (*yfun)(nqy);
      if (!used_q) { 
  *lengthscpy++= run_length; 
  run_length= 1; 
  count++;
      }
    }

    /* remaining points */
    for ( ipt= 2; ipt<npts; ipt++ ) {
      px= qx; py= qy;
      qx= *xarray++; qy= *yarray++;
      if ( clip_one_line( px, py, qx, qy, xmin, xmax, ymin, ymax,
        &npx, &npy, &nqx, &nqy, &used_p, &used_q, &idcode ) ) {
  run_length++;
  if (!used_p) { *ixcpy++= (*xfun)(npx); *iycpy++= (*yfun)(npy); }
  *ixcpy++= (*xfun)(nqx); *iycpy++= (*yfun)(nqy);
  if (!used_q) { 
    *lengthscpy++= run_length;
    count++;
    run_length= 1;
  }
      }
    }

    /* close possibly open run */
    if (used_q) { *lengthscpy= run_length; count++; }
  }
  else {
    (void) no_clip( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
         xfun, yfun );
    count= 1;
    *lengths= npts;
  }

  *pixbuf= ixbuf;
  *piybuf= iybuf;
  *plengths= lengths;
  return( count );
}

static int clean_degenerate_edges(run_length)
int run_length;
/*
This routine scans ixbuf and iybuf, looking for groups of three
points which are colinear in the x or y directions.  The middle
element of such groups is deleted, and the lists are rescanned.
This removes degenerate edges which don't couple distinct non-
degenerate polygons.  Also, any polygon which drops below three
edges is summarily eliminated.  The final number of edges is 
returned.
*/
{
  int old_run_length, changed= 1, i, j;

  old_run_length= run_length;

  while (changed) {
    changed= 0;

    /* Check for not enough vertices */
    if (run_length <= 2) run_length= 0;

    /* Check for colinear triples in x direction */
    for (i=0; i<run_length-2; i++)
      if ( (ixbuf[i]==ixbuf[i+1]) && (ixbuf[i+1]==ixbuf[i+2]) ) {
  for (j=i+1; j<run_length-1; j++) {
    ixbuf[j]= ixbuf[j+1];
    iybuf[j]= iybuf[j+1];
  }
  run_length--;
  changed= 1;
      }

    /* Check for colinear triples in y direction */
    for (i=0; i<run_length-2; i++)
      if ( (iybuf[i]==iybuf[i+1]) && (iybuf[i+1]==iybuf[i+2]) ) {
  for (j=i+1; j<run_length-1; j++) {
    ixbuf[j]= ixbuf[j+1];
    iybuf[j]= iybuf[j+1];
  }
  run_length--;
  changed= 1;
      }
    }

  if (debug) fprintf(stderr,"clean_degenerate_edges:removed %d vertices\n",
         old_run_length - run_length);

  return( run_length );
}

static int add_one_edge( hit, idcode, px, py, qx, qy, npx, npy, nqx, nqy, 
      xmin, xmax, ymin, ymax, xfun, yfun, ixcpy, iycpy )
int hit, idcode, *ixcpy, *iycpy;
float px, py, qx, qy, npx, npy, nqx, nqy, xmin, xmax, ymin, ymax;
int (*xfun)(), (*yfun)();
/*
This routine takes a single set of vertex coordinates and, based on
its idcode (which contains information on where it started and ended-
see clip_one_line()) and whether or not it hit the clipping window,
adds appropriate vertices to the index list to implement the
Liang-Barsky Polygon Clipping Algorithm (see clip_polygon for
reference).  The total number of vertices added is returned.
*/
{
  int count= 0;

  /* The following macros are for use only within this routine */
#define add( x, y ) { *ixcpy++= (*xfun)(x); *iycpy++= (*yfun)(y); count++; }
#define addP add( npx, npy )
#define addQ add( nqx, nqy )
#define addTL add( xmin, ymax )
#define addTR add( xmax, ymax )
#define addBL add( xmin, ymin )
#define addBR add( xmax, ymin )

  switch (idcode) {

  /* p in 1 */
  case 0x99: /* q in 1 */ break;
  case 0x98: /* q in 2 */ break;
  case 0x9a: /* q in 3 */ addTR; break;
  case 0x91: /* q in 4 */ break;
  case 0x90: /* q in 5 */ addP; addQ; break;
  case 0x92: /* q in 6 */ if (hit) { addP; addQ; } else addTR; break;
  case 0x95: /* q in 7 */ addBL; break;
  case 0x94: /* q in 8 */ if (hit) { addP; addQ; } else addBL; break;
  case 0x96: /* q in 9 */ 
    if (hit) { addP; addQ; addBR; }
    else { 
      if leftside( px, py, xmax, ymax, qx, qy ) addTR else addBL;
      addBR;
    }
    break;

  /* p in 2 */
  case 0x89: /* q in 1 */ addTL; break;
  case 0x88: /* q in 2 */ break;
  case 0x8a: /* q in 3 */ addTR; break;
  case 0x81: /* q in 4 */ if (hit) { addP; addQ; } else addTL; break;
  case 0x80: /* q in 5 */ addP; addQ; break;
  case 0x82: /* q in 6 */ if (hit) { addP; addQ; } else addTR; break;
  case 0x85: /* q in 7 */ if (hit) { addP; addQ; } else addTL; addBL; break;
  case 0x84: /* q in 8 */ addP; addQ; break;
  case 0x86: /* q in 9 */ if (hit) { addP; addQ; } else addTR; addBR; break;

  /* p in 3 */
  case 0xa9: /* q in 1 */ addTL; break;
  case 0xa8: /* q in 2 */ break;
  case 0xaa: /* q in 3 */ break;
  case 0xa1: /* q in 4 */ if (hit) { addP; addQ; } else addTL; break;
  case 0xa0: /* q in 5 */ addP; addQ; break;
  case 0xa2: /* q in 6 */ break;
  case 0xa5: /* q in 7 */ 
    if (hit) { addP; addQ; addBL }
    else {
      if leftside( px, py, xmax, ymin, qx, qy ) addBR else addTL;
      addBL;
    }
    break;
  case 0xa4: /* q in 8 */ if (hit) { addP; addQ; } else addBR; break;
  case 0xa6: /* q in 9 */ addBR; break;

  /* p in 4 */
  case 0x19: /* q in 1 */ addTL; break;
  case 0x18: /* q in 2 */ if (hit) { addP; addQ; } else addTL; break;
  case 0x1a: /* q in 3 */ if (hit) { addP; addQ; } else addTL; addTR; break;
  case 0x11: /* q in 4 */ break;
  case 0x10: /* q in 5 */ addP; addQ; break;
  case 0x12: /* q in 6 */ addP; addQ; break;
  case 0x15: /* q in 7 */ addBL; break;
  case 0x14: /* q in 8 */ if (hit) { addP; addQ; } else addBL; break;
  case 0x16: /* q in 9 */ if (hit) { addP; addQ; } else addBL; addBR; break;

  /* p in 5 */
  case 0x09: /* q in 1 */ addQ; addTL; break;
  case 0x08: /* q in 2 */ addQ; break;
  case 0x0a: /* q in 3 */ addQ; addTR; break;
  case 0x01: /* q in 4 */ addQ; break;
  case 0x00: /* q in 5 - case of trivial acceptance */ addQ; break;
  case 0x02: /* q in 6 */ addQ; break;
  case 0x05: /* q in 7 */ addQ; addBL; break;
  case 0x04: /* q in 8 */ addQ; break;
  case 0x06: /* q in 9 */ addQ; addBR; break;

  /* p in 6 */
  case 0x29: /* q in 1 */ if (hit) { addP; addQ; } else addTR; addTL; break;
  case 0x28: /* q in 2 */ if (hit) { addP; addQ; } else addTR; break;
  case 0x2a: /* q in 3 */ addTR; break;
  case 0x21: /* q in 4 */ addP; addQ; break;
  case 0x20: /* q in 5 */ addP; addQ; break;
  case 0x22: /* q in 6 */ break;
  case 0x25: /* q in 7 */ if (hit) { addP; addQ; } else addBR; addBL; break;
  case 0x24: /* q in 8 */ if (hit) { addP; addQ; } else addBR; break;
  case 0x26: /* q in 9 */ addBR; break;

  /* p in 7 */
  case 0x59: /* q in 1 */ addTL; break;
  case 0x58: /* q in 2 */ if (hit) { addP; addQ; } else addTL; break;
  case 0x5a: /* q in 3 */ 
    if (hit) { addP; addQ; addTR }
    else {
      if leftside( px, py, xmin, ymax, qx, qy ) addTL else addBR;
      addTR;
    }
    break;
  case 0x51: /* q in 4 */ break;
  case 0x50: /* q in 5 */ addP; addQ; break;
  case 0x52: /* q in 6 */ if (hit) { addP; addQ; } else addTR; break;
  case 0x55: /* q in 7 */ break;
  case 0x54: /* q in 8 */ break;
  case 0x56: /* q in 9 */ addBR; break;

  /* p in 8 */
  case 0x49: /* q in 1 */ if (hit) { addP; addQ; } else addBL; addTL; break;
  case 0x48: /* q in 2 */ addP; addQ; break;
  case 0x4a: /* q in 3 */ if (hit) { addP; addQ; } else addBR; addTR; break;
  case 0x41: /* q in 4 */ if (hit) { addP; addQ; } else addBL; break;
  case 0x40: /* q in 5 */ addP; addQ; break;
  case 0x42: /* q in 6 */ if (hit) { addP; addQ; } else addBR; break;
  case 0x45: /* q in 7 */ addBL; break;
  case 0x44: /* q in 8 */ break;
  case 0x46: /* q in 9 */ addBR; break;

  /* p in 9 */
  case 0x69: /* q in 1 */ 
    if (hit) { addP; addQ; addTL }
    else {
      if leftside( px, py, xmin, ymin, qx, qy ) addBL else addTR;
      addTL;
    }
    break;
  case 0x68: /* q in 2 */ if (hit) { addP; addQ; } else addTR; break;
  case 0x6a: /* q in 3 */ addTR; break;
  case 0x61: /* q in 4 */ if (hit) { addP; addQ; } else addBL; break;
  case 0x60: /* q in 5 */ addP; addQ; break;
  case 0x62: /* q in 6 */ break;
  case 0x65: /* q in 7 */ addBL; break;
  case 0x64: /* q in 8 */ break;
  case 0x66: /* q in 9 */ break;

  default:
    fprintf(stderr,"invalid code 0x%x in clip_one_line!\n",idcode);
    exit(2);
  }

  return( count );  /* value accumulated by add macro */

  /* Undo macro definition */
#undef addP
#undef addQ
#undef addTL
#undef addTR
#undef addBL
#undef addBR
#undef add
  /* NOTREACHED */
}

int clip_polygon( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
      xfun, yfun, pixbuf, piybuf, plengths )
int npts, **pixbuf, **piybuf, **plengths, (*xfun)(), (*yfun)();
float *xarray, *yarray, xmin, xmax, ymin, ymax;
/*
This routine converts a coordinate list to device (integer) coordinates,
clipping as for a polygon.  The value returned is the total number of
polygons to be plotted;  the number of vertices in each polygon is
stored in the array pointed to by *plengths.

At the moment, only one polygon is returned for each clipped polygon.
This polygon may have degenerate edges, tracing around the edge of the
clipping window to join two separate non-degenerate regions.

The algorithm is basically the Liang-Barsky Polygon Algorithm (see
Computer Graphics Principals and Practice by Foley, van Dam, Feiner,
and Hughes).  This algorithm may produce degenerate edges;  they
will ultimately be cleaned up by the routine 'clean_degenerate_edges'.
Degenerate edges connecting disjoint non-degenerate regions 
of the clipped polygon will survive, however.
*/
{
  int ipt, *ixcpy, *iycpy, *lengthscpy, count, hit;
  int used_p, used_q, idcode, run_length;
  float *xcpy, *ycpy, px, py, qx, qy, npx, npy, nqx, nqy;
  
  if (debug) fprintf(stderr, "clip_polygon:\n");

  if (npts < 3) {
    fprintf(stderr," clip_polygon: got only %d vertices, need at least 3!\n",
      npts);
    return( 0 );
  }

  if ( getcmem( 2*npts, 2*npts, 1 ) != 1 ) {
    fprintf(stderr," Error allocating memory in WRTPGN; call ignored.\n");
    return( 0 );
  };

  
  if (clip_flag) {
    ixcpy= ixbuf; 
    iycpy= iybuf;
    lengthscpy= lengths;
    run_length= 0;
    xcpy= xarray;
    ycpy= yarray;

    /* Go through all the vertices */
    qx= *xcpy++; 
    qy= *ycpy++;
    for ( ipt= 1; ipt<npts; ipt++ ) {
      px= qx; py= qy;
      qx= *xcpy++; qy= *ycpy++;
      hit= clip_one_line( px, py, qx, qy, xmin, xmax, ymin, ymax,
       &npx, &npy, &nqx, &nqy, &used_p, &used_q, &idcode );
      run_length += 
  add_one_edge( hit, idcode, px, py, qx, qy, npx, npy, nqx, nqy, 
         xmin, xmax, ymin, ymax, xfun, yfun, ixcpy, iycpy );
      ixcpy= ixbuf + run_length;
      iycpy= iybuf + run_length;
    }

    /* Pick up the connection between the last vertex and the first */
    px= qx; py= qy;
    qx= *xarray; qy= *yarray;
    hit= clip_one_line( px, py, qx, qy, xmin, xmax, ymin, ymax,
           &npx, &npy, &nqx, &nqy, &used_p, &used_q, &idcode );
    run_length += 
      add_one_edge( hit, idcode, px, py, qx, qy, npx, npy, nqx, nqy, 
       xmin, xmax, ymin, ymax, xfun, yfun, ixcpy, iycpy );

  /* The following routine removes degenerate edges.  It needs no
   * parameters, as it works on ixbuf, iybuf, and lengths directly.
   */
  run_length= clean_degenerate_edges(run_length);

  count= 1;
  *lengths= run_length;
  }
  else {
    (void) no_clip( npts, xarray, yarray, xmin, xmax, ymin, ymax, 
         xfun, yfun );
    count= 1;
    *lengths= npts;
  }

  *pixbuf= ixbuf;
  *piybuf= iybuf;
  *plengths= lengths;
  return(count);
}