polygon_intersection_gpc.cpp

这是一个GPS相关的程序

    /* Redirect any p->proxy references to q->proxy */
    for (target= p->proxy; list; list= list->next)
    {
      if (list->proxy == target)
      {
        list->active= FALSE;
        list->proxy= q->proxy;
      }
    }
  }
}


static void add_local_min(polygon_node **p, edge_node *edge,
                          double x, double y)
{
  polygon_node *existing_min;
  vertex_node  *nv;

  existing_min= *p;

  _MALLOC(*p, sizeof(polygon_node), "polygon node creation");

  /* Create a new vertex node and set its fields */
  _MALLOC(nv, sizeof(vertex_node), "vertex node creation");
  nv->x= x;
  nv->y= y;
  nv->next= NULL;

  /* Initialise proxy to point to p itself */
  (*p)->proxy= (*p);
  (*p)->active= TRUE;
  (*p)->next= existing_min;

  /* Make v[LEFT] and v[RIGHT] point to new vertex nv */
  (*p)->v[LEFT]= nv;
  (*p)->v[RIGHT]= nv;

  /* Assign polygon p to the edge */
  edge->outp[ABOVE]= *p;
}


static int count_tristrips(polygon_node *tn)
{
  int total;

  for (total= 0; tn; tn= tn->next)
    if (tn->active > 2)
      total++;
  return total;
}


static void add_vertex(vertex_node **t, double x, double y)
{
  if (!(*t))
  {
    _MALLOC(*t, sizeof(vertex_node), "tristrip vertex creation");
    (*t)->x= x;
    (*t)->y= y;
    (*t)->next= NULL;
  }
  else
    /* Head further down the list */
    add_vertex(&((*t)->next), x, y);
}


static void new_tristrip(polygon_node **tn, edge_node *edge,
                         double x, double y)
{
  if (!(*tn))
  {
    _MALLOC(*tn, sizeof(polygon_node), "tristrip node creation");
    (*tn)->next= NULL;
    (*tn)->v[LEFT]= NULL;
    (*tn)->v[RIGHT]= NULL;
    (*tn)->active= 1;
    add_vertex(&((*tn)->v[LEFT]), x, y); 
    edge->outp[ABOVE]= *tn;
  }
  else
    /* Head further down the list */
    new_tristrip(&((*tn)->next), edge, x, y);
}


static bbox *create_contour_bboxes(gpc_polygon *p)
{
  bbox *box;
  int   c, v;

  _MALLOC(box, p->num_contours * sizeof(bbox), "Bounding box creation");

  /* Construct contour bounding boxes */
  for (c= 0; c < p->num_contours; c++)
  {
    /* Initialise bounding box extent */
    box[c].xmin= DBL_MAX;
    box[c].ymin= DBL_MAX;
    box[c].xmax= -DBL_MAX;
    box[c].ymax= -DBL_MAX;

    for (v= 0; v < p->contour[c].num_vertices; v++)
    {
      /* Adjust bounding box */
      if (p->contour[c].vertex[v].x < box[c].xmin)
        box[c].xmin= p->contour[c].vertex[v].x;
      if (p->contour[c].vertex[v].y < box[c].ymin)
        box[c].ymin= p->contour[c].vertex[v].y;
      if (p->contour[c].vertex[v].x > box[c].xmax)
        box[c].xmax= p->contour[c].vertex[v].x;
      if (p->contour[c].vertex[v].y > box[c].ymax)
          box[c].ymax= p->contour[c].vertex[v].y;
    }
  }
  return box;  
}


static void minimax_test(gpc_polygon *subj, gpc_polygon *clip, gpc_op op)
{
  bbox *s_bbox, *c_bbox;
  int   s, c, *o_table, overlap;

  s_bbox= create_contour_bboxes(subj);
  c_bbox= create_contour_bboxes(clip);

  _MALLOC(o_table, subj->num_contours * clip->num_contours * sizeof(int),
         "overlap table creation");

  /* Check all subject contour bounding boxes against clip boxes */
  for (s= 0; s < subj->num_contours; s++)
    for (c= 0; c < clip->num_contours; c++)
      o_table[c * subj->num_contours + s]=
             (!((s_bbox[s].xmax < c_bbox[c].xmin) ||
                (s_bbox[s].xmin > c_bbox[c].xmax))) &&
             (!((s_bbox[s].ymax < c_bbox[c].ymin) ||
                (s_bbox[s].ymin > c_bbox[c].ymax)));

  /* For each clip contour, search for any subject contour overlaps */
  for (c= 0; c < clip->num_contours; c++)
  {
    overlap= 0;
    for (s= 0; (!overlap) && (s < subj->num_contours); s++)
      overlap= o_table[c * subj->num_contours + s];

    if (!overlap)
      /* Flag non contributing status by negating vertex count */
      clip->contour[c].num_vertices = -clip->contour[c].num_vertices;
  }  

  if (op == GPC_INT)
  {  
    /* For each subject contour, search for any clip contour overlaps */
    for (s= 0; s < subj->num_contours; s++)
    {
      overlap= 0;
      for (c= 0; (!overlap) && (c < clip->num_contours); c++)
        overlap= o_table[c * subj->num_contours + s];

      if (!overlap)
        /* Flag non contributing status by negating vertex count */
        subj->contour[s].num_vertices = -subj->contour[s].num_vertices;
    }  
  }

  FREE(s_bbox);
  FREE(c_bbox);
  FREE(o_table);
}


/*
===========================================================================
                             Public Functions
===========================================================================
*/

void gpc_free_polygon(gpc_polygon *p)
{
  int c;

  for (c= 0; c < p->num_contours; c++)
    FREE(p->contour[c].vertex);
  FREE(p->hole);
  FREE(p->contour);
  p->num_contours= 0;
}


void gpc_read_polygon(FILE *fp, int read_hole_flags, gpc_polygon *p)
{
  int c, v;

  fscanf(fp, "%d", &(p->num_contours));
  _MALLOC(p->hole, p->num_contours * sizeof(int),
         "hole flag array creation");
  _MALLOC(p->contour, p->num_contours
         * sizeof(gpc_vertex_list), "contour creation");
  for (c= 0; c < p->num_contours; c++)
  {
    fscanf(fp, "%d", &(p->contour[c].num_vertices));

    if (read_hole_flags)
      fscanf(fp, "%d", &(p->hole[c]));
    else
      p->hole[c]= FALSE; /* Assume all contours to be external */

    _MALLOC(p->contour[c].vertex, p->contour[c].num_vertices
           * sizeof(gpc_vertex), "vertex creation");
    for (v= 0; v < p->contour[c].num_vertices; v++)
      fscanf(fp, "%lf %lf", &(p->contour[c].vertex[v].x),
                            &(p->contour[c].vertex[v].y));
  }
}


void gpc_write_polygon(FILE *fp, int write_hole_flags, gpc_polygon *p)
{
  int c, v;

  fprintf(fp, "%d\n", p->num_contours);
  for (c= 0; c < p->num_contours; c++)
  {
    fprintf(fp, "%d\n", p->contour[c].num_vertices);

    if (write_hole_flags)
      fprintf(fp, "%d\n", p->hole[c]);
    
    for (v= 0; v < p->contour[c].num_vertices; v++)
      fprintf(fp, "% .*lf % .*lf\n",
              DBL_DIG, p->contour[c].vertex[v].x,
              DBL_DIG, p->contour[c].vertex[v].y);
  }
}


void gpc_add_contour(gpc_polygon *p, gpc_vertex_list *new_contour, int hole)
{
  int             *extended_hole, c, v;
  gpc_vertex_list *extended_contour;

  /* Create an extended hole array */
  _MALLOC(extended_hole, (p->num_contours + 1)
         * sizeof(int), "contour hole addition");

  /* Create an extended contour array */
  _MALLOC(extended_contour, (p->num_contours + 1)
         * sizeof(gpc_vertex_list), "contour addition");

  /* Copy the old contour and hole data into the extended arrays */
  for (c= 0; c < p->num_contours; c++)
  {
    extended_hole[c]= p->hole[c];
    extended_contour[c]= p->contour[c];
  }

  /* Copy the new contour and hole onto the end of the extended arrays */
  c= p->num_contours;
  extended_hole[c]= hole;
  extended_contour[c].num_vertices= new_contour->num_vertices;
  _MALLOC(extended_contour[c].vertex, new_contour->num_vertices
         * sizeof(gpc_vertex), "contour addition");
  for (v= 0; v < new_contour->num_vertices; v++)
    extended_contour[c].vertex[v]= new_contour->vertex[v];

  /* Dispose of the old contour */
  FREE(p->contour);
  FREE(p->hole);

  /* Update the polygon information */
  p->num_contours++;
  p->hole= extended_hole;
  p->contour= extended_contour;
}


void gpc_polygon_clip(gpc_op op, gpc_polygon *subj, gpc_polygon *clip,
                      gpc_polygon *result)
{
  sb_tree       *sbtree= NULL;
  it_node       *it= NULL, *intersect;
  edge_node     *edge, *prev_edge, *next_edge, *succ_edge, *e0, *e1;
  edge_node     *aet= NULL, *c_heap= NULL, *s_heap= NULL;
  lmt_node      *lmt= NULL, *local_min;
  polygon_node  *out_poly= NULL, *p, *q, *poly, *npoly, *cf= NULL;
  vertex_node   *vtx, *nv;
  h_state        horiz[2];
  int            in[2], exists[2], parity[2]= {LEFT, LEFT};
  int            c, v, contributing, search, scanbeam= 0, sbt_entries= 0;
  int            vclass, bl, br, tl, tr;
  double        *sbt= NULL, xb, px, yb, yt, dy, ix, iy;

  /* Test for trivial NULL result cases */
  if (((subj->num_contours == 0) && (clip->num_contours == 0))
   || ((subj->num_contours == 0) && ((op == GPC_INT) || (op == GPC_DIFF)))
   || ((clip->num_contours == 0) &&  (op == GPC_INT)))
  {
    result->num_contours= 0;
    result->hole= NULL;
    result->contour= NULL;
    return;
  }

  /* Identify potentialy contributing contours */
  if (((op == GPC_INT) || (op == GPC_DIFF))
   && (subj->num_contours > 0) && (clip->num_contours > 0))
    minimax_test(subj, clip, op);

  /* Build LMT */
  if (subj->num_contours > 0)
    s_heap= build_lmt(&lmt, &sbtree, &sbt_entries, subj, SUBJ, op);
  if (clip->num_contours > 0)
    c_heap= build_lmt(&lmt, &sbtree, &sbt_entries, clip, CLIP, op);

  /* Return a NULL result if no contours contribute */
  if (lmt == NULL)
  {
    result->num_contours= 0;
    result->hole= NULL;
    result->contour= NULL;
    reset_lmt(&lmt);
    FREE(s_heap);
    FREE(c_heap);
    return;
  }

  /* Build scanbeam table from scanbeam tree */
  _MALLOC(sbt, sbt_entries * sizeof(double), "sbt creation");
  build_sbt(&scanbeam, sbt, sbtree);
  scanbeam= 0;
  free_sbtree(&sbtree);

  /* Allow pointer re-use without causing memory leak */
  if (subj == result)
    gpc_free_polygon(subj);
  if (clip == result)
    gpc_free_polygon(clip);

  /* Invert clip polygon for difference operation */
  if (op == GPC_DIFF)
    parity[CLIP]= RIGHT;

  local_min= lmt;

  /* Process each scanbeam */
  while (scanbeam < sbt_entries)
  {
    /* Set yb and yt to the bottom and top of the scanbeam */
    yb= sbt[scanbeam++];
    if (scanbeam < sbt_entries)
    {
      yt= sbt[scanbeam];
      dy= yt - yb;
    }

    /* === SCANBEAM BOUNDARY PROCESSING ================================ */

    /* If LMT node corresponding to yb exists */
    if (local_min)
    {
      if (local_min->y == yb)
      {
        /* Add edges starting at this local minimum to the AET */
        for (edge= local_min->first_bound; edge; edge= edge->next_bound)
          add_edge_to_aet(&aet, edge, NULL);

        local_min= local_min->next;
      }
    }

    /* Set dummy previous x value */
    px= -DBL_MAX;

    /* Create bundles within AET */
    e0= aet;
    e1= aet;

    /* Set up bundle fields of first edge */
    aet->bundle[ABOVE][ aet->type]= (aet->top.y != yb);
    aet->bundle[ABOVE][!aet->type]= FALSE;
    aet->bstate[ABOVE]= UNBUNDLED;

    for (next_edge= aet->next; next_edge; next_edge= next_edge->next)
    {
      /* Set up bundle fields of next edge */
      next_edge->bundle[ABOVE][ next_edge->type]= (next_edge->top.y != yb);
      next_edge->bundle[ABOVE][!next_edge->type]= FALSE;
      next_edge->bstate[ABOVE]= UNBUNDLED;

      /* Bundle edges above the scanbeam boundary if they coincide */
      if (next_edge->bundle[ABOVE][next_edge->type])
      {
        if (EQ(e0->xb, next_edge->xb) && EQ(e0->dx, next_edge->dx)
	 && (e0->top.y != yb))
        {
          next_edge->bundle[ABOVE][ next_edge->type]^= 
            e0->bundle[ABOVE][ next_edge->type];
          next_edge->bundle[ABOVE][!next_edge->type]= 
            e0->bundle[ABOVE][!next_edge->type];
          next_edge->bstate[ABOVE]= BUNDLE_HEAD;
          e0->bundle[ABOVE][CLIP]= FALSE;
          e0->bundle[ABOVE][SUBJ]= FALSE;
          e0->bstate[ABOVE]= BUNDLE_TAIL;
        }
        e0= next_edge;
      }
    }
    
    horiz[CLIP]= NH;
    horiz[SUBJ]= NH;

    /* Process each edge at this scanbeam boundary */
    for (edge= aet; edge; edge= edge->next)
    {
      exists[CLIP]= edge->bundle[ABOVE][CLIP] + 
                   (edge->bundle[BELOW][CLIP] << 1);
      exists[SUBJ]= edge->bundle[ABOVE][SUBJ] + 
                   (edge->bundle[BELOW][SUBJ] << 1);