#define         max(x,y)  ((x) > (y) ? (x) : (y))
#define         min(x,y)  ((x) < (y) ? (x) : (y))

#define         ERROR   0.00000001      /* ZERO approximation *

#define INSIDE 1
#define OUTSIDE 2
#define INTERSECT 3

/* position of point with respect to polygon. retruns INSIDE, OUTSIDE or INTERSECT
cannot handle situations such that described in "Algorithms", R. Sedgewick, pp. 355 */

/*
points: number of polygon's vertices
point:  array which holds the coordinates of the polygon's vertices
p:      the point we try to determine if it is inside, outside or on 
        the polygon.

struct POINT {
float x,y;
};

struct LINE_SEGMENT {
   struct POINT *p1, *p2;
};

*/

int point_to_polygon(p, point, points)
POINT *p, *point;
int points;
{
        int i, active;
        POINT r, p1, p2, *pptr;
        LINE_SEGMENT horiz, line_seg;
        float x1, x2, y1, y2;
        float active_xs[MAXPOINTS];
 
        p1.x = 0;  p1.y = p->y;
        p2.x = MAXNUMBER;  p2.y = p->y;
        horiz.p1 = &p1; horiz.p2 = &p2;
        active = 0;
 
        pptr = point;
        for (i = 0; i < points - 1; i++) {
          y1 = min(pptr->y, (pptr + 1)->y);
          y2 = max(pptr->y, (pptr + 1)->y);
          if (y1 < p->y && p->y <= y2) {
            line_seg.p1 = pptr;
            line_seg.p2 = pptr + 1;
            intersection(&horiz, &line_seg, &r); /* finds intersection point */
            active++;
            active_xs[active - 1] = r.x;
          }
          pptr++;
        }
          
        if (active > 1) {
          sort(active_xs, active); /* sorts active_xs, active: array size */
          for (i = 0; i < active - 1 ; i++) if (i%2 == 0) {
            x1 = min(active_xs[i], active_xs[i + 1]);
            x2 = max(active_xs[i], active_xs[i + 1]);
            if (ZERO(x1 - p->x) || ZERO(x2 - p->x)) return(INTERSECT);
            if ( (x1 < p->x) && (p->x < x2) ) return(INSIDE);
          }
        }   
        return(OUTSIDE);
}


/* returns in p the intersection of two line segments. if parallel returns 0 else
returns 1 */

int intersection(line_seg1, line_seg2, p)
LINE_SEGMENT *line_seg1, *line_seg2;
POINT *p;
{
        register float x1,y1, x2,y2, x3,y3, x4,y4;
        register float A1, B1, C1, A2, B2, C2;
        
        x1 = (*line_seg1).p1->x;  y1 = (*line_seg1).p1->y;
        x2 = (*line_seg1).p2->x;  y2 = (*line_seg1).p2->y;

        x3 = (*line_seg2).p1->x;  y3 = (*line_seg2).p1->y;
        x4 = (*line_seg2).p2->x;  y4 = (*line_seg2).p2->y;

        if (ZERO(x1 - x2)) {
           A1 = 1; B1 = 0; C1 = -x1;
        } else {
           A1 = (y2 - y1)/(x2 - x1); B1 = -1; C1 = y1 - A1*x1;
        }
 
        if  (ZERO(x3 - x4)) {
          A2 = 1; B2 = 0; C2 = -x3;
        } else {
          A2 = (y4 - y3)/(x4 - x3); B2 = -1; C2 = y3 - A2*x3;
        }
          
        if (ZERO(A1*B2 - B1*A2)) return(0); /* lines parallel */
        else {
          (*p).x = (-C1*B2 + C2*B1)/(A1*B2 - A2*B1);
          (*p).y = (-C2*A1 + C1*A2)/(A1*B2 - A2*B1);
          return(1)
        }
}  

/* ZERO approximation */
int ZERO(number)
float number;
{
        if (fabs((double)number) < ERROR) return(1);
        else return(0);
}