home *** CD-ROM | disk | FTP | other *** search
/ Current Shareware 1994 January / SHAR194.ISO / graphuti / polyops.zip / POLYCALC.C < prev    next >
Text File  |  1991-07-19  |  27KB  |  871 lines
  1.  
  2.  
  3. /* Partial implementation in C of the program described in "An algorithm
  4.    for computing the union,intersection and difference of two polygons",
  5.    A. Margalit and G. D. Knott, Comput. & Graphics, 13, 2, 167-183, 1989.
  6.  
  7.    Only the case of simple, regular curves is handled.
  8.  
  9.    Programmed by:
  10.    B. J. Ball, 3304 Glen Rose Drive, Austin, Texas 78731 (CompuServe 71141,2321)
  11. */
  12.  
  13.  
  14. #include <stdio.h>
  15. #include <math.h>
  16.  
  17. #define MAXPTS   100    /* maximum size of ANY polygon used, not just input */
  18. #define MAXX     1.1    /* greater than any permissible x-coordinate */
  19. #define MAXCOMPS  20    /* maximum number of output polygons allowed */
  20.  
  21. #define MIN(a,b)  (a) < (b) ? (a) : (b)
  22. #define MAX(a,b)  (a) > (b) ? (a) : (b)
  23.  
  24. typedef struct {double x; double y;} point;
  25. typedef struct {point first; point second;} segment;
  26. typedef struct {int nverts; double x[MAXPTS]; double y[MAXPTS];} polygon;
  27. typedef struct {point v; int class; int next;} ventry; /* for vertex rings */
  28. typedef struct {point p1; point p2; int del;} fentry;  /* for frags list */
  29.  
  30. /* function prototypes */
  31. void
  32. changeorientation(polygon *),
  33. consolidate(int *,double x[],double y[]),
  34. initializeAV_BV(void),
  35. insertAV(point,int,int),
  36. insertBV(point,int,int),
  37. insertE(point,point);
  38.  
  39. int
  40. build_poly(int),
  41. cross_check(point *,segment,segment),
  42. equal(point,point),
  43. equalflt(double,double),
  44. findintersection(segment *,segment,segment),
  45. findposA(point,segment,int),
  46. findposB(point,segment,int),
  47. gt(double,double),
  48. gte(double,double),
  49. inside(polygon,point),
  50. onseg(point,segment,double *),
  51. orientation(polygon),
  52. polygonsoperation(void),
  53. simple(polygon);
  54.  
  55. /* prompt() is part of the demo program, used here to display error messages;
  56.    it should be replaced or rewritten if polycalc.c is used independently. */
  57. extern int prompt(int, int, int, char *);
  58.  
  59. ventry AV[MAXPTS],BV[MAXPTS];      /* linked lists of vertices and isects */
  60. fentry EF[MAXPTS];                 /* edge fragment list */
  61. polygon outpoly[MAXCOMPS];         /* program output */
  62. int ptype[MAXCOMPS];               /* types of the output polygons */
  63. int numcomps;                      /* number of output polygons */
  64. int skip[MAXPTS];                  /* flags used in build_poly() function */
  65.  
  66. polygon A,B;                       /* input polygons */
  67. int oper;                /* operation: 0,1,2,3 = intersection,union,A-B,B-A */
  68. int Atype=0, Btype=0;              /* 0 = "island", 1 = "hole" */
  69.  
  70. int AVsize, BVsize, EFsize;        /* current sizes of AV,BV,EF arrays */
  71.  
  72. double roundoff = .0000001;   /* for checking equality of float calculations */
  73.  
  74. /* Specialized arrays used for distinguishing various cases */
  75. int polygonsorientation[2][2][4] =
  76.             {{{1,1,-1,-1},{-1,-1,1,1}},{{-1,-1,1,1},{1,1,-1,-1}}};
  77. int fragmenttype[2][2][4][2] =
  78.           {{{{1,1},{-1,-1},{-1,1},{1,-1}},{{-1,-1},{1,-1},{1,1},{-1,-1}}},
  79.            {{{1,-1},{-1,1},{-1,-1},{1,1}},{{-1,-1},{1,1},{1,-1},{-1,1}}}};
  80. int resultorientation[2][2][4] =
  81.           {{{1,1,1,-1},{1,-1,1,1}},{{-1,1,1,1},{1,1,-1,1}}};
  82.  
  83.  
  84. /* ====================================================================== */
  85. /*                  principal routine                                     */
  86. /* ====================================================================== */
  87. /*
  88.   Input:  polygons A,B
  89.           int oper      (operation: 0=intersection, 1=union, 2=A-B, 3=B-A)
  90.   Output: int numcomp = number of output components
  91.           polygon outpoly[numcomps] = array of output polygons
  92.           int ptype[numcomps] = array giving type of each output polygon
  93. */
  94.  
  95. polygonsoperation(void)
  96. {
  97.    int orientationA, orientationB;
  98.    int i,i1,j,j1,k,type,t1,t2,cnt,k1,k2,err_flag=0;
  99.    double u1,u2;
  100.    segment a,b,c,tmp;
  101.    point m,p1,p2;
  102.  
  103.    orientationA = orientation(A);
  104.    orientationB = orientation(B);
  105.    if (polygonsorientation[Atype][Btype][oper] == 1)
  106.      {
  107.       if (orientationA != orientationB)
  108.          changeorientation(&B);
  109.      }
  110.    else if (orientationA == orientationB)
  111.       changeorientation(&B);
  112.    initializeAV_BV();
  113.  
  114.  /* find intersections of sides of A with sides of B, insert in VA and VB */
  115.    for (i=0; i<A.nverts; i++)
  116.      {
  117.       i1 = (i+1) % A.nverts;
  118.       a.first.x = AV[i].v.x, a.first.y = AV[i].v.y;
  119.       a.second.x = AV[i1].v.x, a.second.y = AV[i1].v.y;
  120.       for (j=0; j<B.nverts; j++)
  121.         {
  122.          j1 = (j+1) % B.nverts;
  123.          b.first.x = BV[j].v.x, b.first.y = BV[j].v.y;
  124.          b.second.x = BV[j1].v.x, b.second.y = BV[j1].v.y;
  125.          k = findintersection(&c,a,b);
  126.          if (k==1)                      /* only one point of intersection */
  127.            {
  128.             k1 = findposA(c.first,a,i); /* last pt. of AV[] <= c.first on a */
  129.             insertAV(c.first,0,k1);     /* insert c.first here */
  130.             k1 = findposB(c.first,b,j); /* repeat for BV */
  131.             insertBV(c.first,0,k1);
  132.            }
  133.           else if (k==2)                 /* two points of intersection */
  134.            {
  135.             /* set p1,p2 to c.first,c.second in order on a */
  136.             onseg(c.first,a,&u1);   /* u1,u2 are coordinates of.. */
  137.             onseg(c.second,a,&u2);  /* .. c.first,c.second on segment a */
  138.             if (gt(u2,u1))
  139.               {
  140.                p1.x = c.first.x, p1.y = c.first.y;
  141.                p2.x = c.second.x, p2.y = c.second.y;
  142.               }
  143.             else
  144.               {
  145.                p2.x = c.first.x, p2.y = c.first.y;
  146.                p1.x = c.second.x, p1.y = c.second.y;
  147.               }
  148.             k1 = findposA(p1,a,i);         /* last point of AV[] <= p1 on a */
  149.             if (equal(AV[k1].v,p1))        /* p1,p2 should follow in order */
  150.               insertAV(p2,0,k1);
  151.             else
  152.               {
  153.                insertAV(p1,0,k1);           /* set AV[AVsize] = p1, AVsize++ */
  154.                insertAV(p2,0,AVsize-1);     /* next position following p1 */
  155.               }
  156.         /* repeat above procedure for B */
  157.             onseg(c.first,b,&u1);   /* u1,u2 are coordinates of.. */
  158.             onseg(c.second,b,&u2);  /* .. c.first,c.second on segment b */
  159.            if (gt(u2,u1))
  160.               {
  161.                p1.x = c.first.x, p1.y = c.first.y;
  162.                p2.x = c.second.x, p2.y = c.second.y;
  163.               }
  164.             else
  165.               {
  166.                p2.x = c.first.x, p2.y = c.first.y;
  167.                p1.x = c.second.x, p1.y = c.second.y;
  168.               }
  169.             k1 = findposB(p1,b,j);         /* last point of BV[] <= p1 on a */
  170.             if (equal(BV[k1].v,p1))        /* p1,p2 should follow in order */
  171.               insertBV(p2,0,k1);
  172.             else
  173.               {
  174.                insertBV(p1,0,k1);           /* set BV[AVsize] = p1, BVsize++ */
  175.                insertBV(p2,0,BVsize-1);     /* next position following p1 */
  176.               }
  177.            }
  178.         }
  179.      }
  180.  
  181. /* select edge fragments and insert into EF[] */
  182.    type = fragmenttype[Atype][Btype][oper][0];   /* last index specifies A */
  183.  
  184.    i = 0; j = AV[i].next; cnt=0;
  185.    while (cnt < AVsize)
  186.      {
  187.       t1 = AV[i].class;
  188.       t2 = AV[j].class;
  189.       if (t1 == type || t2 == type)
  190.          insertE(AV[i].v,AV[j].v);
  191.       else if (t1==0 && t2==0)
  192.         {
  193.          m.x = (AV[i].v.x + AV[j].v.x)/2;
  194.          m.y = (AV[i].v.y + AV[j].v.y)/2;
  195.          k = inside(B,m);
  196.  
  197.          if (k==type || k==0)
  198.             insertE(AV[i].v,AV[j].v);
  199.         }
  200.       i = j; j = AV[i].next; cnt++;
  201.      }
  202.  
  203.    type = fragmenttype[Atype][Btype][oper][1];   /* last index specifies B */
  204.  
  205.    i = 0; j = BV[i].next; cnt = 0;
  206.    while (cnt < BVsize)
  207.      {
  208.       t1 = BV[i].class;
  209.       t2 = BV[j].class;
  210.       if (t1 == type || t2 == type)
  211.          insertE(BV[i].v,BV[j].v);
  212.       if (t1==0 && t2==0)
  213.         {
  214.          m.x = (BV[i].v.x + BV[j].v.x)/2;
  215.          m.y = (BV[i].v.y + BV[j].v.y)/2;
  216.          k = inside(A,m);
  217.          if (k==type || k==0)
  218.             insertE(BV[i].v,BV[j].v);
  219.         }
  220.       i = j; j = BV[i].next; cnt++;
  221.      }
  222.  
  223. /* arrange edge fragments into polygons */
  224.    k = 0;
  225.    while ((err_flag = build_poly(k)) > 0)
  226.      {
  227.       i = orientation(outpoly[k]);
  228.       if (i == orientationA)
  229.         {
  230.          if (resultorientation[Atype][Btype][oper] == 1)
  231.             ptype[k] = Atype;
  232.          else
  233.             ptype[k] = 1-Atype;
  234.         }
  235.       else
  236.         {
  237.          if (resultorientation[Atype][Btype][oper] == 1)
  238.             ptype[k] = 1-Atype;
  239.          else
  240.             ptype[k] = Atype;
  241.         }
  242.        k++;
  243.        if (k > MAXCOMPS)
  244.          {
  245.           err_flag = -1;
  246.           prompt(32,192,15,
  247.             "Too many components in result. Press any key to continue");
  248.           getch();
  249.           break;
  250.          }
  251.      }
  252.    numcomps = k;          /* number of output polygons found */
  253.    if (err_flag)
  254.       return -1;          /* build_poly failed */
  255.    else
  256.       return numcomps;
  257. }
  258.  
  259. /* ------------------ initializeAV_BV() ----------------------------------- */
  260. /* classify original vertices and insert them in vertex rings AV[],BV[]     */
  261. void
  262. initializeAV_BV(void)
  263. {
  264.    int i,n;
  265.    point p;
  266.  
  267.    n = A.nverts;
  268.    for (i=0; i<n; i++)
  269.      {
  270.       p.x = AV[i].v.x = A.x[i];
  271.       p.y = AV[i].v.y = A.y[i];
  272.       AV[i].class = inside(B,p);
  273.       AV[i].next = i+1;
  274.      }
  275.    AV[n-1].next = 0.0;         /* close up ring */
  276.    AVsize = n;
  277.  
  278.    n = B.nverts;
  279.    for (i=0; i<n; i++)
  280.      {
  281.       p.x = BV[i].v.x = B.x[i];
  282.       p.y = BV[i].v.y = B.y[i];
  283.       BV[i].class = inside(A,p);
  284.       BV[i].next = i+1;
  285.       }
  286.    BV[n-1].next = 0.0;         /* close up ring */
  287.    BVsize = n;
  288. }
  289.  
  290.  
  291. /* ======================================================================== */
  292. /*      Principal supporting routines, listed alphabetically                */
  293. /* ======================================================================== */
  294.  
  295. /* --------------- construct a polygon from boundary fragments ------------- */
  296. /* fills outpoly[k], marks points used */
  297.  
  298. int
  299. build_poly(int k)
  300. {
  301.    int i,j,j1,nverts,match_found;
  302.    point u0,u1;   /* u0 = first vertex of outpoly, u1 = current last vertex */
  303.  
  304.    for (i=0; i<EFsize; i++)
  305.      {
  306.       if (skip[i] || EF[i].del)  continue;
  307.       skip[i] = 1;
  308.       outpoly[k].x[0] = u0.x = EF[i].p1.x;
  309.       outpoly[k].y[0] = u0.y = EF[i].p1.y;  /* u0 = EF[i].p1 */
  310.       outpoly[k].x[1] = u1.x = EF[i].p2.x;
  311.       outpoly[k].y[1] = u1.y = EF[i].p2.y;  /* u1 = EF[i].p2 */
  312.       nverts = 2;       /* number of vertices added to  outpoly[k] so far */
  313.       do
  314.         {
  315.          match_found = 0;                         /* continue flag */
  316.          j1 = 0;                                  /* index of matching seg */
  317.          for (j=i+1; j<EFsize; j++)
  318.            {
  319.             if (skip[j] || EF[j].del)  continue;
  320.             if (equal(EF[j].p1,u1))               /* if EF[j].p1 = u1, then.. */
  321.               {
  322.                j1 = j;
  323.                if (equal(EF[j].p2,u0))       /* if back at starting point ... */
  324.                  {
  325.                   skip[j] = 1;
  326.                   outpoly[k].nverts = nverts;    /* .. polygon completed */
  327.                              /* eliminate collinear sequences of vertices */
  328.                   consolidate(&outpoly[k].nverts,outpoly[k].x,outpoly[k].y);
  329.                   return 1;
  330.                  }
  331.               }
  332.            }
  333.          if (!j1)
  334.            {
  335.             prompt(32,192,15,
  336.                  "Computation failed - press any key to continue");
  337.             getch();
  338.             goto err_exit;
  339.            }
  340.          skip[j1] = 1;
  341.          u1.x = EF[j1].p2.x, u1.y = EF[j1].p2.y;   /* ..set u1 = EF[j].p2 */
  342.          match_found = 1;               /* continue flag */
  343.          outpoly[k].x[nverts] = u1.x;
  344.          outpoly[k].y[nverts] = u1.y;
  345.          nverts++;
  346.         } while (nverts <= 2*MAXPTS && match_found);
  347.  
  348.       if (nverts > 2*MAXPTS)
  349.         {
  350.          prompt(32,192,15,
  351.          "Too many sides in one component. Press any key to continue");
  352.          getch();
  353.          goto err_exit;
  354.         }
  355.      }
  356.    return 0;
  357. err_exit:
  358.    return -1;
  359. }
  360.  
  361. /* ----------------  change orientation ------------------------------- */
  362. /* (reverse order of vertices of polygon) */
  363. void
  364. changeorientation(polygon *pgn)
  365. {
  366.    int i,j,n = pgn->nverts;
  367.    double tempx[MAXPTS], tempy[MAXPTS];
  368.  
  369.    for (i=0; i<n; i++)
  370.      {
  371.       tempx[i] = pgn->x[i];
  372.       tempy[i] = pgn->y[i];
  373.      }
  374.    for (i=0,j=n-1; i<n; i++,j--)
  375.      {
  376.       pgn->x[i] = tempx[j];
  377.       pgn->y[i] = tempy[j];
  378.      }
  379. }
  380.  
  381. /* --------------------------------------------------------------------- */
  382. /* collinear() returns 1 if the input points are collinear, 0 otherwise  */
  383. int
  384. collinear(double a1, double a2, double b1, double b2, double c1, double c2)
  385. {
  386.    if (fabs((a1-b1)*(a2-c2) - (a1-c1)*(a2-b2)) > roundoff)
  387.       return 0;
  388.    else
  389.       return 1;
  390. }
  391.  
  392. /* ----------------------------------------------------------------------- */
  393. /* routine for eliminating collinear sequences of points from x[0],...x[n] */
  394. /* resets n to new value if necessary */
  395. void
  396. consolidate(int *n, double x[], double y[])
  397. {
  398.    int i=0,j,m=*n;
  399.  
  400.    while (i < m-2)
  401.      {
  402.       if (!collinear(x[i],y[i],x[i+1],y[i+1], x[i+2],y[i+2]))
  403.          i++;
  404.       else
  405.         {
  406.          for (j=i+1; j<m-1; j++)     /* delete middle point, (x[i+1],y[i+1]) */
  407.             x[j]=x[j+1], y[j]=y[j+1];
  408.          m--;
  409.         }
  410.      }
  411.    if (m > 2 && collinear(x[m-2],y[m-2],x[m-1],y[m-1],x[0],y[0]))
  412.       m--;                            /* delete last point, if necessary */
  413.    *n = m;          /* change original value of n to the calculated value */
  414. }
  415.  
  416. /* ---------------- cross_check(point *,segment,segment) ------------------- */
  417. /* returns -1 if segments a and b are parallel, 0 if they do not intersect,
  418.    and 1 if the intersection is a single point. In this last case, the point
  419.    of intersection is returned in p. This routine is from "A Programmer's
  420.    Geometry", Bowyer and Woodwark, London, 1988. */
  421. int
  422. cross_check(point *p, segment a, segment b)
  423. {
  424.    double xlk = a.second.x - a.first.x, ylk = a.second.y - a.first.y;
  425.    double xnm = b.second.x - b.first.x, ynm = b.second.y - b.first.y;
  426.    double xmk = b.first.x - a.first.x, ymk = b.first.y - a.first.y;
  427.    double det,detinv,s,t,x,y;
  428.  
  429.    det = xnm*ylk - ynm*xlk;
  430.    if (fabs(det) < roundoff)
  431.       return -1;
  432.    detinv = 1.0/det;
  433.    s = (xnm*ymk - ynm*xmk)*detinv;
  434.    if (gt(0,s) || gt(s,1))
  435.      return 0;
  436.    t = (xlk*ymk - ylk*xmk)*detinv;
  437.    if (gt(0,t) || gt(t,1))
  438.      return 0;
  439.    p->x = a.first.x + xlk*s;
  440.    p->y = a.first.y + ylk*s;
  441.    return 1;
  442. }
  443.  
  444. /* ---------------- find intersection of two segments ----------------------- */
  445. /* Input: segment a, segment b, segment *c.
  446.    Returns 0 if the intersection of a and b is empty, 1 if it is a single point
  447.    and 2 if it is a nondegenerate segment. In the latter two cases, the segment
  448.    c is set to the intersection of a and b. */
  449.  
  450. int
  451. findintersection(segment *c, segment a, segment b)
  452. {
  453.    int i,cnt;
  454.    point p[2];
  455.    double tmp;          /* dummy variable for onseg() routine */
  456.  
  457.    /* run simple box tests first */
  458.    if (gt(MIN(b.first.x,b.second.x),MAX(a.first.x,a.second.x)))
  459.      return 0;
  460.    if (gt(MIN(a.first.x,a.second.x),MAX(b.first.x,b.second.x)))
  461.      return 0;
  462.    if (gt(MIN(b.first.y,b.second.y),MAX(a.first.y,a.second.y)))
  463.      return 0;
  464.    if (gt(MIN(a.first.y,a.second.y),MAX(b.first.y,b.second.y)))
  465.      return 0;
  466.    /* if box tests do not rule out intersection, use cross_check() routine */
  467.    i = cross_check(&p[0],a,b);  /* test for crossing segments */
  468.    if (i==0)                    /* no intersection */
  469.       return 0;
  470.    if (i==1)                     /* only one point of intersection */
  471.      {
  472.       c->first.x = c->second.x = p[0].x;
  473.       c->first.y = c->second.y = p[0].y;
  474.       return 1;
  475.      }
  476.    cnt = 0;         /* check endpoints, quit when two intersections found */
  477.    if (onseg(a.first,b,&tmp))
  478.      {
  479.       p[cnt].x = a.first.x;
  480.       p[cnt].y = a.first.y;
  481.       cnt++;
  482.      }
  483.    if (onseg(a.second,b,&tmp))
  484.      {
  485.       p[cnt].x = a.second.x;
  486.       p[cnt].y = a.second.y;
  487.       cnt++;
  488.       if (cnt>1)
  489.          goto b;
  490.      }
  491.    if (onseg(b.first,a,&tmp))
  492.      {
  493.       p[cnt].x = b.first.x;
  494.       p[cnt].y = b.first.y;
  495.       cnt++;
  496.       if (cnt>1)
  497.          goto b;
  498.      }
  499.    if (onseg(b.second,a,&tmp))
  500.      {
  501.       p[cnt].x = b.second.x;
  502.       p[cnt].y = b.second.y;
  503.       cnt++;
  504.       if (cnt>1)
  505.          goto b;
  506.      }
  507.  
  508.    if (!cnt)
  509.       return 0;
  510.    else
  511.      {
  512.       c->first.x = c->second.x = p[0].x;
  513.       c->first.y = c->second.y = p[0].y;
  514.       return 1;
  515.      }
  516.  
  517. b: c->first.x = p[0].x, c->first.y = p[0].y;
  518.    c->second.x = p[1].x, c->second.y = p[1].y;
  519.    return 2;
  520. }
  521.  
  522. /* ---------------- find order of added vertices on polygon sides --------- */
  523. /* Input: point p, segment a containing p, integer j; for findposA, a = AV[j]
  524.    and for findposB(), p = BV[j]. Return the index of the last point in AV[]
  525.    or BV[], respectively, which does not follow p on segment a. */
  526.  
  527. int
  528. findposA(point p, segment a, int j)
  529. {
  530.    int i,i0;
  531.    double t,t0,tmp;
  532.  
  533.    if (!onseg(p,a,&t))
  534.      {
  535.       textcolor(15);
  536.       cputs("internal error - cannot continue");
  537.       getch();
  538.       return -1;
  539.      }
  540.    i0 = j; t0 = 0.0;
  541.    for (i=0; i<AVsize; i++)
  542.      {
  543.       if (!onseg(AV[i].v,a,&tmp))  continue;
  544.       if (equalflt(tmp,t))
  545.         {
  546.          i0=i;
  547.          break;
  548.         }
  549.       if (gt(tmp,t0) && gt(t,tmp))      /* t0 < tmp < t */
  550.          t0 = tmp, i0 = i;
  551.      }
  552.    return i0;
  553. }
  554.  
  555. /* ---------------------------------- */
  556. int
  557. findposB(point p, segment a, int j)
  558. {
  559.    int i,i0;
  560.    double t,t0,tmp;
  561.  
  562.    if (!onseg(p,a,&t))
  563.      {
  564.       textcolor(15);
  565.       cputs("internal error - cannot continue");
  566.       getch();
  567.       return -1;
  568.      }
  569.  
  570.    i0 = j; t0 = 0.0;
  571.    for (i=0; i<BVsize; i++)
  572.      {
  573.       if (!onseg(BV[i].v,a,&tmp))  continue;
  574.       if (equalflt(tmp,t))
  575.         {
  576.          i0=i;
  577.          break;
  578.         }
  579.       if (gt(tmp,t0) && gt(t,tmp))      /* t0 < tmp < t */
  580.          t0 = tmp, i0 = i;
  581.      }
  582.    return i0;
  583. }
  584.  
  585. /* ---------------- insertAV(point,type,k) -------------------------------- */
  586. /* add point at end of AV, then change links so new point follows current k */
  587. void
  588. insertAV(point p, int type, int k)
  589. {
  590.    int i,n;
  591.  
  592.    /* check on current appearances of p in AV[] */
  593.    for (i=0; i<AVsize; i++)
  594.       if (equal(AV[i].v, p))
  595.          return;         /* a point can appear at most once */
  596.    n = AVsize;
  597.    AVsize++;
  598.    AV[n].v.x = p.x, AV[n].v.y = p.y;
  599.    AV[n].next = AV[k].next;
  600.    AV[n].class = type;
  601.    AV[k].next = n;
  602. }
  603.  
  604. /* ---------------- insertBV(point,type,k) -------------------------------- */
  605. /* add point at end of BV, then change links so new point follows current k */
  606. void
  607. insertBV(point p, int type, int k)
  608. {
  609.    int i,n;
  610.  
  611.    /* check on current appearances of p in BV[] */
  612.    for (i=0; i<BVsize; i++)
  613.       if (equal(BV[i].v, p))
  614.          return;         /* a point can appear at most once */
  615.    n = BVsize;
  616.    BVsize++;
  617.    BV[n].v.x = p.x, BV[n].v.y = p.y;
  618.    BV[n].next = BV[k].next;
  619.    BV[n].class = type;
  620.    BV[k].next = n;
  621. }
  622.  
  623. /* ---------------- insertE(point,point) ---------------------------------- */
  624. /* add segment to list of edge fragments if it is not already there, except */
  625. /* if the inverse of the segment in in EF, delete it and don't insert seg   */
  626. /* (if non-regular output is to be allowed, this routine must be modified)  */
  627. void
  628. insertE(point a, point b)
  629. {
  630.    int i,k1,k2;
  631.  
  632. /* check for presence of segment ab or ba in EF[] */
  633.    for (i=0; i<EFsize; i++)
  634.     {
  635.      k1 = ((equal(EF[i].p1,a) && equal(EF[i].p2,b)));
  636.      if (k1)                         /* ab found in EF[] */
  637.         return;
  638.      k2 = (equal(EF[i].p1,b) && equal(EF[i].p2,a));
  639.      if (k2)                          /* ba found in EF[] */
  640.        {
  641.         EF[i].del = 1;
  642.         return;
  643.        }
  644.     }
  645.    EF[EFsize].p1.x = a.x, EF[EFsize].p1.y = a.y;
  646.    EF[EFsize].p2.x = b.x, EF[EFsize].p2.y = b.y;
  647.    EFsize++;
  648. }
  649.  
  650. /* ---------------------- inside(polygon,point) -------------------------- */
  651. /* returns -1,0,1 according as point p is outside, on boundary, inside pgn */
  652. int
  653. inside(polygon pgn, point p)
  654. {
  655.    int i,j,k,cnt;
  656.    double y1,tmp;    /* tmp is dummy variable for onseg() routine */
  657.    point p1,p2;
  658.    segment a,b,c;
  659.  
  660.    /* set segment a to start at p and go horizontally right beyond pgn */
  661.    a.first.x = p.x, a.first.y = a.second.y = p.y; a.second.x = MAXX;
  662.  
  663.    /* check intersection of segment a with each side of pgn */
  664.    cnt = 0;
  665.    for (i=0; i<pgn.nverts; i++)
  666.      {
  667.       b.first.x = pgn.x[i], b.first.y = pgn.y[i];
  668.       j = (i+1) % pgn.nverts;
  669.       b.second.x = pgn.x[j], b.second.y = pgn.y[j];
  670.       if (onseg(p,b,&tmp))
  671.          return 0;
  672.       if (gt(p.x,MAX(b.first.x,b.second.x))
  673.           || gt(p.y,MAX(b.first.y,b.second.y))
  674.              || gt(MIN(b.first.y,b.second.y),p.y))
  675.         continue;                              /* b does not intersect a */
  676.       k = findintersection(&c,a,b);
  677.       if (k)
  678.         {
  679.          y1 = MIN(b.first.y, b.second.y);
  680.          if (!equalflt(c.first.y,y1) && !equalflt(c.second.y,y1))
  681.             cnt++;
  682.         }
  683.       }
  684.     if (cnt %2)
  685.        return 1;
  686.     else
  687.        return -1;
  688. }
  689.  
  690. /* ----------------------- onseg(p,s,&t) ---------------------------------- */
  691. /* returns 0 if point is not on segment; otherwise returns 1 and sets t to
  692.    the coordinate of p on s; i.e., p = (s.first) + t*(s.second)             */
  693. int
  694. onseg(point p, segment seg, double *t)
  695. {
  696.    double x1=seg.first.x, y1=seg.first.y, x2=seg.second.x, y2=seg.second.y;
  697.    double x=p.x, y=p.y;
  698.    double temp;
  699.  
  700.    if (equalflt(x1,x2))                  /* vertical segment */
  701.      {
  702.       if (!equalflt(x,x1))
  703.          return 0;
  704.       if (equalflt(y1,y2))               /* degenerate segment */
  705.         {
  706.          if (equalflt(y,y1)) return 1;
  707.          else return 0;
  708.         }
  709.       temp = (y-y1)/(y2-y1);
  710.       goto a;
  711.      }
  712.    if (equalflt(y1,y2))                  /* horizontal segment */
  713.      {
  714.       if (!equalflt(y,y1))
  715.          return 0;
  716.       if (equalflt(x1,x2))               /* degenerate segment */
  717.         {
  718.          if (equalflt(x,x1)) return 1;
  719.          else return 0;
  720.         }
  721.       temp = (x-x1)/(x2-x1);
  722.       goto a;
  723.      }
  724.    if (!equalflt((x-x1)*(y2-y1),(y-y1)*(x2-x1)))
  725.       return 0;
  726.    temp = (x-x1)/(x2-x1);
  727. a:
  728.    *t = temp;
  729.    if (gte(temp,0.0) && gte(1.0,temp))    /* if 0 <= temp <= 1, p is on seg */
  730.       return 1;
  731.    else
  732.       return 0;
  733. }
  734.  
  735. /* ---------------------- orientation(polygon) ----------------------------- */
  736. /*      returns +1 for counterclockwise orientation, -1 for clockwise        */
  737. /* (the routine described in the Margalit-Knott paper does not work in all   */
  738. /* situations which the present program presents to it; there is a kludge    */
  739. /* here which makes things all right in the current program but still fails  */
  740. /* for non-regular polygons. Clearly the authors had something in mind other */
  741. /* than what is given in the paper, but it is not clear exactly what the     */
  742. /* routine ought to be in the general case. Surely not this ugly mess.)      */
  743.  
  744. int
  745. orientation(polygon pgn)
  746. {
  747.    point p,q,r;
  748.    int i,j,k,l,m,n=pgn.nverts,i0,i1,j0=0,temp;
  749.  
  750.    if (n<=2) return 0;      /* empty or degenerate polygon */
  751.  
  752.    /* find a vertex of pgn with minimum x-coordinate */
  753.    for (j=0; j<n; j++)
  754.       if (pgn.x[j] < pgn.x[j0])
  755.          j0 = j;
  756.    q.x = pgn.x[j0], q.y = pgn.y[j0];   /* q is a left-most vertex of pgn */
  757.    if (j0>0)
  758.       i0 = j0-1;
  759.    else
  760.       i0 = n-1;   /* i0 is the index of an immediate predecessor of q on pgn */
  761.    for (i=0; i<n-1; i++)
  762.      {
  763.       i1 = (i+1)%n;
  764.       if (pgn.x[i1]==q.x && pgn.y[i+1]==q.y && pgn.x[i]<pgn.x[i0])
  765.          i0 = i;
  766.      }
  767.    i = i0, j = (i+1)%n, k = (j+1) % n;
  768.    /* now i is the index of the left-most immediate predecessor of q,
  769.           j is an index of q in pgn, and
  770.           k is the index of the next point on pgn */
  771.  
  772.    p.x = pgn.x[i], p.y = pgn.y[i];
  773.    r.x = pgn.x[k], r.y = pgn.y[k];
  774.  
  775.  /* if either of the segments pq, qr is vertical, the orientation is obvious */
  776.     if (equalflt(p.x, q.x))
  777.       if (gt(p.y,q.y))
  778.          return 1;
  779.       else
  780.          return -1;
  781.  
  782.    if (equalflt(q.x, r.x))
  783.       if (gt(q.y,r.y))
  784.          return 1;
  785.       else
  786.          return -1;
  787. /* otherwise, calculate the orientation from the slopes of pq and qr */
  788.    if (gt((p.y-q.y)/(p.x-q.x), (q.y-r.y)/(q.x-r.x)))
  789.       temp = 1;
  790.    else
  791.       temp = -1;
  792.  
  793. /* if q is an interior point of a  side other than pq and qp, that side must
  794.  be vertical; if this situation obtains, the calculation just made is wrong.
  795.  Check it out. */
  796.  /* since this routine was written, the build_poly() function has been changed
  797.     so that the situation envisioned here should not occur. However ... */
  798.    for (l=0; l<n; l++)
  799.      {
  800.       if (l==i || l==j || l==k)  continue;
  801.       if (!equalflt(pgn.x[l], q.x)) continue;
  802.       m = (l+1) % n;
  803.       if (!equalflt(pgn.x[m], q.x)) continue;
  804.       if ((gt(q.y,pgn.y[l]) && gt(pgn.y[m],q.y))
  805.            || (gt(q.y,pgn.y[m]) && gt(pgn.y[l],q.y)))
  806.               return -temp;
  807.      }
  808.    return temp;
  809. }
  810.  
  811. /* ----------------------------------------------------------------------- */
  812. /* check input polygons for simplicity */
  813. int
  814. simple(polygon pgn)
  815. {
  816.    int i,i1,j,j1,n=pgn.nverts;
  817.    segment c, seg[MAXPTS];
  818.  
  819.    for (i=0; i<n; i++)
  820.      {
  821.       i1 = (i+1)%n;
  822.       seg[i].first.x = pgn.x[i];
  823.       seg[i].first.y = pgn.y[i];
  824.       seg[i].second.x = pgn.x[i1];
  825.       seg[i].second.y = pgn.y[i1];
  826.      }
  827.  
  828.    for (j=2; j<n-1; j++)
  829.      if (findintersection(&c,seg[0],seg[j]))
  830.         return 0;
  831.    for (i=1; i<n-2; i++)
  832.       for (j=i+2; j<n; j++)
  833.         if (findintersection(&c,seg[i],seg[j]))
  834.            return 0;
  835.   return 1;
  836. }
  837.  
  838. /* ======================================================================= */
  839. /*              Miscellaneous floating point utilities                     */
  840. /* ======================================================================= */
  841. /* the value of 'roundoff' might have to be adjusted for some applications */
  842.  
  843.       /* check two doubles for equality */
  844. int
  845. equalflt(double x, double y)
  846. {
  847.    return (fabs(x-y)<roundoff);
  848. }
  849.  
  850.       /* check for x > y */
  851. int
  852. gt(double x, double y)
  853. {
  854.    return (x-y > roundoff);
  855. }
  856.  
  857.      /* check for x >= y */
  858. int
  859. gte(double x, double y)
  860. {
  861.    return (gt(x,y) + equalflt(x,y));
  862. }
  863.  
  864.      /* check two points for equality */
  865. int
  866. equal(point a, point b)
  867. {
  868.    return ((fabs(a.x-b.x)<roundoff) && (fabs(a.y-b.y)<roundoff));
  869. }
  870.  
  871.