Ticket #754: polygon_hole.diff

File polygon_hole.diff, 99.8 KB (added by tegzed, 9 years ago)

fixed intersection detection for horizontal and vertical lines

  • maptool/maptool.h

     
    121121struct boundary {
    122122        struct item_bin *ib;
    123123        struct country_table *country;
     124        enum item_type* area_item_types;
    124125        char *iso2;
    125126        GList *segments,*sorted_segments;
    126127        GList *children;
     
    187188GList *geom_poly_segments_sort(GList *in, enum geom_poly_segment_type type);
    188189struct geom_poly_segment *item_bin_to_poly_segment(struct item_bin *ib, int type);
    189190int geom_poly_segments_point_inside(GList *in, struct coord *c);
     191GList* geom_poly_segments_group(GList *in, GList *out);
    190192void clip_line(struct item_bin *ib, struct rect *r, struct tile_parameter *param, struct item_bin_sink *out);
    191193void clip_polygon(struct item_bin *ib, struct rect *r, struct tile_parameter *param, struct item_bin_sink *out);
     194int self_intersect_test(struct item_bin*ib);
     195int geom_point_in_polygon(struct coord*c, struct geom_poly_segment*poly);
     196int geom_polygon_in_polygon(struct geom_poly_segment*poly1, struct geom_poly_segment*poly2);
     197int geom_polygon_is_clockwise(struct geom_poly_segment*poly);
     198void geom_polygon_reverse(struct geom_poly_segment*poly);
    192199
    193200/* itembin.c */
    194201
     
    240247extern struct buffer node_buffer;
    241248extern int processed_nodes, processed_nodes_out, processed_ways, processed_relations, processed_tiles;
    242249extern struct item_bin *item_bin;
     250extern struct item_bin *item_bin_relation_area;
    243251extern int bytes_read;
    244252extern int overlap;
    245253extern int unknown_country;
  • maptool/itembin_buffer.c

     
    2424
    2525static char buffer[800000];
    2626struct item_bin *item_bin=(struct item_bin *)(void *)buffer;
     27static char buffer_relation_area[8000000];
     28struct item_bin *item_bin_relation_area=(struct item_bin *)(void *)buffer_relation_area;
    2729static struct node_item *node_item=(struct node_item *)(void *)buffer;
    2830
    2931struct node_item *
  • maptool/seidel/tri.c

     
     1#include "triangulate.h"
     2#include <sys/time.h>
     3
     4
     5static int initialise(n)
     6     int n;
     7{
     8  register int i;
     9
     10  for (i = 1; i <= n; i++)
     11    seg[i].is_inserted = FALSE;
     12
     13  generate_random_ordering(n);
     14 
     15  return 0;
     16}
     17
     18#ifdef STANDALONE
     19
     20int main(argc, argv)
     21     int argc;
     22     char *argv[];
     23{
     24  int n, nmonpoly, genus;
     25  int op[SEGSIZE][3], i, ntriangles;
     26
     27  if ((argc < 2) || ((n = read_segments(argv[1], &genus)) < 0))
     28    {
     29      fprintf(stderr, "usage: triangulate <filename>\n");
     30      exit(1);
     31    }
     32
     33  initialise(n);
     34  construct_trapezoids(n);
     35  nmonpoly = monotonate_trapezoids(n);
     36  ntriangles = triangulate_monotone_polygons(n, nmonpoly, op);
     37
     38  for (i = 0; i < ntriangles; i++)
     39    printf("triangle #%d: %d %d %d\n", i,
     40           op[i][0], op[i][1], op[i][2]);
     41
     42  return 0;
     43}
     44
     45
     46#else  /* Not standalone. Use this as an interface routine */
     47
     48
     49/* Input specified as contours.
     50 * Outer contour must be anti-clockwise.
     51 * All inner contours must be clockwise.
     52 * 
     53 * Every contour is specified by giving all its points in order. No
     54 * point shoud be repeated. i.e. if the outer contour is a square,
     55 * only the four distinct endpoints shopudl be specified in order.
     56 * 
     57 * ncontours: #contours
     58 * cntr: An array describing the number of points in each
     59 *       contour. Thus, cntr[i] = #points in the i'th contour.
     60 * vertices: Input array of vertices. Vertices for each contour
     61 *           immediately follow those for previous one. Array location
     62 *           vertices[0] must NOT be used (i.e. i/p starts from
     63 *           vertices[1] instead. The output triangles are
     64 *           specified  w.r.t. the indices of these vertices.
     65 * triangles: Output array to hold triangles.
     66 * 
     67 * Enough space must be allocated for all the arrays before calling
     68 * this routine
     69 */
     70
     71
     72int triangulate_polygon(ncontours, cntr, vertices, triangles)
     73     s32 ncontours;
     74     s32 cntr[];
     75     s32 (*vertices)[2];
     76     s32 (*triangles)[3];
     77{
     78  register int i;
     79  int nmonpoly, ccount, npoints, genus;
     80  int n;
     81
     82  memset((void *)seg, 0, sizeof(seg));
     83  ccount = 0;
     84  i = 1;
     85  while (ccount < ncontours)
     86    {
     87      int j;
     88      int first, last;
     89
     90      npoints = cntr[ccount];
     91      first = i;
     92      last = first + npoints - 1;
     93      for (j = 0; j < npoints; j++, i++)
     94        {
     95          seg[i].v0.x = (double)(vertices[i][0]);
     96          seg[i].v0.y = (double)vertices[i][1];
     97
     98          if (i == last)
     99            {
     100              seg[i].next = first;
     101              seg[i].prev = i-1;
     102              seg[i-1].v1 = seg[i].v0;
     103            }
     104          else if (i == first)
     105            {
     106              seg[i].next = i+1;
     107              seg[i].prev = last;
     108              seg[last].v1 = seg[i].v0;
     109            }
     110          else
     111            {
     112              seg[i].prev = i-1;
     113              seg[i].next = i+1;
     114              seg[i-1].v1 = seg[i].v0;
     115            }
     116         
     117          seg[i].is_inserted = FALSE;
     118        }
     119     
     120      ccount++;
     121    }
     122 
     123  genus = ncontours - 1;
     124  n = i-1;
     125
     126  initialise(n);
     127  if(construct_trapezoids(n)) {
     128    return 0;
     129  }
     130  nmonpoly = monotonate_trapezoids(n);
     131  return triangulate_monotone_polygons(n, nmonpoly, triangles);
     132}
     133
     134
     135/* This function returns TRUE or FALSE depending upon whether the
     136 * vertex is inside the polygon or not. The polygon must already have
     137 * been triangulated before this routine is called.
     138 * This routine will always detect all the points belonging to the
     139 * set (polygon-area - polygon-boundary). The return value for points
     140 * on the boundary is not consistent!!!
     141 */
     142
     143int is_point_inside_polygon(vertex)
     144     double vertex[2];
     145{
     146  point_t v;
     147  int trnum, rseg;
     148  trap_t *t;
     149
     150  v.x = vertex[0];
     151  v.y = vertex[1];
     152 
     153  trnum = locate_endpoint(&v, &v, 1);
     154  t = &tr[trnum];
     155 
     156  if (t->state == ST_INVALID)
     157    return FALSE;
     158 
     159  if ((t->lseg <= 0) || (t->rseg <= 0))
     160    return FALSE;
     161  rseg = t->rseg;
     162  return _greater_than_equal_to(&seg[rseg].v1, &seg[rseg].v0);
     163}
     164
     165
     166#endif /* STANDALONE */
  • maptool/seidel/monotone.c

     
     1#include "triangulate.h"
     2#include <math.h>
     3
     4#define CROSS_SINE(v0, v1) ((v0).x * (v1).y - (v1).x * (v0).y)
     5#define LENGTH(v0) (sqrt((v0).x * (v0).x + (v0).y * (v0).y))
     6
     7
     8static monchain_t mchain[TRSIZE]; /* Table to hold all the monotone */
     9                                  /* polygons . Each monotone polygon */
     10                                  /* is a circularly linked list */
     11
     12static vertexchain_t vert[SEGSIZE]; /* chain init. information. This */
     13                                    /* is used to decide which */
     14                                    /* monotone polygon to split if */
     15                                    /* there are several other */
     16                                    /* polygons touching at the same */
     17                                    /* vertex  */
     18
     19static int mon[SEGSIZE];        /* contains position of any vertex in */
     20                                /* the monotone chain for the polygon */
     21static int visited[TRSIZE];
     22static int chain_idx, op_idx, mon_idx;
     23
     24
     25static int triangulate_single_polygon(int, int, int, int (*)[3]);
     26static int traverse_polygon(int, int, int, int);
     27
     28/* Function returns TRUE if the trapezoid lies inside the polygon */
     29static int inside_polygon(t)
     30     trap_t *t;
     31{
     32  int rseg = t->rseg;
     33
     34  if (t->state == ST_INVALID)
     35    return 0;
     36
     37  if ((t->lseg <= 0) || (t->rseg <= 0))
     38    return 0;
     39 
     40  if (((t->u0 <= 0) && (t->u1 <= 0)) ||
     41      ((t->d0 <= 0) && (t->d1 <= 0))) /* triangle */
     42    return (_greater_than(&seg[rseg].v1, &seg[rseg].v0));
     43 
     44  return 0;
     45}
     46
     47
     48/* return a new mon structure from the table */
     49static int newmon()
     50{
     51  return ++mon_idx;
     52}
     53
     54
     55/* return a new chain element from the table */
     56static int new_chain_element()
     57{
     58  return ++chain_idx;
     59}
     60
     61
     62static double get_angle(vp0, vpnext, vp1)
     63     point_t *vp0;
     64     point_t *vpnext;
     65     point_t *vp1;
     66{
     67  point_t v0, v1;
     68 
     69  v0.x = vpnext->x - vp0->x;
     70  v0.y = vpnext->y - vp0->y;
     71
     72  v1.x = vp1->x - vp0->x;
     73  v1.y = vp1->y - vp0->y;
     74
     75  if (CROSS_SINE(v0, v1) >= 0)  /* sine is positive */
     76    return DOT(v0, v1)/LENGTH(v0)/LENGTH(v1);
     77  else
     78    return (-1.0 * DOT(v0, v1)/LENGTH(v0)/LENGTH(v1) - 2);
     79}
     80
     81
     82/* (v0, v1) is the new diagonal to be added to the polygon. Find which */
     83/* chain to use and return the positions of v0 and v1 in p and q */
     84static int get_vertex_positions(v0, v1, ip, iq)
     85     int v0;
     86     int v1;
     87     int *ip;
     88     int *iq;
     89{
     90  vertexchain_t *vp0, *vp1;
     91  register int i;
     92  double angle, temp;
     93  int tp, tq;
     94
     95  vp0 = &vert[v0];
     96  vp1 = &vert[v1];
     97 
     98  /* p is identified as follows. Scan from (v0, v1) rightwards till */
     99  /* you hit the first segment starting from v0. That chain is the */
     100  /* chain of our interest */
     101 
     102  angle = -4.0;
     103  for (i = 0; i < 4; i++)
     104    {
     105      if (vp0->vnext[i] <= 0)
     106        continue;
     107      if ((temp = get_angle(&vp0->pt, &(vert[vp0->vnext[i]].pt),
     108                            &vp1->pt)) > angle)
     109        {
     110          angle = temp;
     111          tp = i;
     112        }
     113    }
     114
     115  *ip = tp;
     116
     117  /* Do similar actions for q */
     118
     119  angle = -4.0;
     120  for (i = 0; i < 4; i++)
     121    {
     122      if (vp1->vnext[i] <= 0)
     123        continue;     
     124      if ((temp = get_angle(&vp1->pt, &(vert[vp1->vnext[i]].pt),
     125                            &vp0->pt)) > angle)
     126        {
     127          angle = temp;
     128          tq = i;
     129        }
     130    }
     131
     132  *iq = tq;
     133
     134  return 0;
     135}
     136
     137 
     138/* v0 and v1 are specified in anti-clockwise order with respect to
     139 * the current monotone polygon mcur. Split the current polygon into
     140 * two polygons using the diagonal (v0, v1)
     141 */
     142static int make_new_monotone_poly(mcur, v0, v1)
     143     int mcur;
     144     int v0;
     145     int v1;
     146{
     147  int p, q, ip, iq;
     148  int mnew = newmon();
     149  int i, j, nf0, nf1;
     150  vertexchain_t *vp0, *vp1;
     151 
     152  vp0 = &vert[v0];
     153  vp1 = &vert[v1];
     154
     155  get_vertex_positions(v0, v1, &ip, &iq);
     156
     157  p = vp0->vpos[ip];
     158  q = vp1->vpos[iq];
     159
     160  /* At this stage, we have got the positions of v0 and v1 in the */
     161  /* desired chain. Now modify the linked lists */
     162
     163  i = new_chain_element();      /* for the new list */
     164  j = new_chain_element();
     165
     166  mchain[i].vnum = v0;
     167  mchain[j].vnum = v1;
     168
     169  mchain[i].next = mchain[p].next;
     170  mchain[mchain[p].next].prev = i;
     171  mchain[i].prev = j;
     172  mchain[j].next = i;
     173  mchain[j].prev = mchain[q].prev;
     174  mchain[mchain[q].prev].next = j;
     175
     176  mchain[p].next = q;
     177  mchain[q].prev = p;
     178
     179  nf0 = vp0->nextfree;
     180  nf1 = vp1->nextfree;
     181
     182  vp0->vnext[ip] = v1;
     183
     184  vp0->vpos[nf0] = i;
     185  vp0->vnext[nf0] = mchain[mchain[i].next].vnum;
     186  vp1->vpos[nf1] = j;
     187  vp1->vnext[nf1] = v0;
     188
     189  vp0->nextfree++;
     190  vp1->nextfree++;
     191
     192#ifdef DEBUG
     193  fprintf(stderr, "make_poly: mcur = %d, (v0, v1) = (%d, %d)\n",
     194          mcur, v0, v1);
     195  fprintf(stderr, "next posns = (p, q) = (%d, %d)\n", p, q);
     196#endif
     197
     198  mon[mcur] = p;
     199  mon[mnew] = i;
     200  return mnew;
     201}
     202
     203/* Main routine to get monotone polygons from the trapezoidation of
     204 * the polygon.
     205 */
     206
     207int monotonate_trapezoids(n)
     208     int n;
     209{
     210  register int i;
     211  int tr_start;
     212
     213  memset((void *)vert, 0, sizeof(vert));
     214  memset((void *)visited, 0, sizeof(visited));
     215  memset((void *)mchain, 0, sizeof(mchain));
     216  memset((void *)mon, 0, sizeof(mon));
     217 
     218  /* First locate a trapezoid which lies inside the polygon */
     219  /* and which is triangular */
     220  for (i = 0; i < TRSIZE; i++)
     221    if (inside_polygon(&tr[i]))
     222      break;
     223  tr_start = i;
     224 
     225  /* Initialise the mon data-structure and start spanning all the */
     226  /* trapezoids within the polygon */
     227
     228#if 0
     229  for (i = 1; i <= n; i++)
     230    {
     231      mchain[i].prev = i - 1;
     232      mchain[i].next = i + 1;
     233      mchain[i].vnum = i;
     234      vert[i].pt = seg[i].v0;
     235      vert[i].vnext[0] = i + 1; /* next vertex */
     236      vert[i].vpos[0] = i;      /* locn. of next vertex */
     237      vert[i].nextfree = 1;
     238    }
     239  mchain[1].prev = n;
     240  mchain[n].next = 1;
     241  vert[n].vnext[0] = 1;
     242  vert[n].vpos[0] = n;
     243  chain_idx = n;
     244  mon_idx = 0;
     245  mon[0] = 1;                   /* position of any vertex in the first */
     246                                /* chain  */
     247
     248#else
     249
     250  for (i = 1; i <= n; i++)
     251    {
     252      mchain[i].prev = seg[i].prev;
     253      mchain[i].next = seg[i].next;
     254      mchain[i].vnum = i;
     255      vert[i].pt = seg[i].v0;
     256      vert[i].vnext[0] = seg[i].next; /* next vertex */
     257      vert[i].vpos[0] = i;      /* locn. of next vertex */
     258      vert[i].nextfree = 1;
     259    }
     260
     261  chain_idx = n;
     262  mon_idx = 0;
     263  mon[0] = 1;                   /* position of any vertex in the first */
     264                                /* chain  */
     265
     266#endif
     267 
     268  /* traverse the polygon */
     269  if (tr[tr_start].u0 > 0)
     270    traverse_polygon(0, tr_start, tr[tr_start].u0, TR_FROM_UP);
     271  else if (tr[tr_start].d0 > 0)
     272    traverse_polygon(0, tr_start, tr[tr_start].d0, TR_FROM_DN);
     273 
     274  /* return the number of polygons created */
     275  return newmon();
     276}
     277
     278
     279/* recursively visit all the trapezoids */
     280static int traverse_polygon(mcur, trnum, from, dir)
     281     int mcur;
     282     int trnum;
     283     int from;
     284     int dir;
     285{
     286  trap_t *t = &tr[trnum];
     287  int howsplit, mnew;
     288  int v0, v1, v0next, v1next;
     289  int retval, tmp;
     290  int do_switch = FALSE;
     291
     292  if ((trnum <= 0) || visited[trnum])
     293    return 0;
     294
     295  visited[trnum] = TRUE;
     296 
     297  /* We have much more information available here. */
     298  /* rseg: goes upwards   */
     299  /* lseg: goes downwards */
     300
     301  /* Initially assume that dir = TR_FROM_DN (from the left) */
     302  /* Switch v0 and v1 if necessary afterwards */
     303
     304
     305  /* special cases for triangles with cusps at the opposite ends. */
     306  /* take care of this first */
     307  if ((t->u0 <= 0) && (t->u1 <= 0))
     308    {
     309      if ((t->d0 > 0) && (t->d1 > 0)) /* downward opening triangle */
     310        {
     311          v0 = tr[t->d1].lseg;
     312          v1 = t->lseg;
     313          if (from == t->d1)
     314            {
     315              do_switch = TRUE;
     316              mnew = make_new_monotone_poly(mcur, v1, v0);
     317              traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     318              traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);     
     319            }
     320          else
     321            {
     322              mnew = make_new_monotone_poly(mcur, v0, v1);
     323              traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     324              traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
     325            }
     326        }
     327      else
     328        {
     329          retval = SP_NOSPLIT;  /* Just traverse all neighbours */
     330          traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     331          traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     332          traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     333          traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     334        }
     335    }
     336 
     337  else if ((t->d0 <= 0) && (t->d1 <= 0))
     338    {
     339      if ((t->u0 > 0) && (t->u1 > 0)) /* upward opening triangle */
     340        {
     341          v0 = t->rseg;
     342          v1 = tr[t->u0].rseg;
     343          if (from == t->u1)
     344            {
     345              do_switch = TRUE;
     346              mnew = make_new_monotone_poly(mcur, v1, v0);
     347              traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     348              traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);     
     349            }
     350          else
     351            {
     352              mnew = make_new_monotone_poly(mcur, v0, v1);
     353              traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     354              traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     355            }
     356        }
     357      else
     358        {
     359          retval = SP_NOSPLIT;  /* Just traverse all neighbours */
     360          traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     361          traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     362          traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     363          traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     364        }
     365    }
     366 
     367  else if ((t->u0 > 0) && (t->u1 > 0))
     368    {
     369      if ((t->d0 > 0) && (t->d1 > 0)) /* downward + upward cusps */
     370        {
     371          v0 = tr[t->d1].lseg;
     372          v1 = tr[t->u0].rseg;
     373          retval = SP_2UP_2DN;
     374          if (((dir == TR_FROM_DN) && (t->d1 == from)) ||
     375              ((dir == TR_FROM_UP) && (t->u1 == from)))
     376            {
     377              do_switch = TRUE;
     378              mnew = make_new_monotone_poly(mcur, v1, v0);
     379              traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     380              traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     381              traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     382              traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     383            }
     384          else
     385            {
     386              mnew = make_new_monotone_poly(mcur, v0, v1);
     387              traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     388              traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     389              traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     390              traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);       
     391            }
     392        }
     393      else                      /* only downward cusp */
     394        {
     395          if (_equal_to(&t->lo, &seg[t->lseg].v1))
     396            {
     397              v0 = tr[t->u0].rseg;
     398              v1 = seg[t->lseg].next;
     399
     400              retval = SP_2UP_LEFT;
     401              if ((dir == TR_FROM_UP) && (t->u0 == from))
     402                {
     403                  do_switch = TRUE;
     404                  mnew = make_new_monotone_poly(mcur, v1, v0);
     405                  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     406                  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     407                  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     408                  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
     409                }
     410              else
     411                {
     412                  mnew = make_new_monotone_poly(mcur, v0, v1);
     413                  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     414                  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     415                  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     416                  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     417                }
     418            }
     419          else
     420            {
     421              v0 = t->rseg;
     422              v1 = tr[t->u0].rseg;     
     423              retval = SP_2UP_RIGHT;
     424              if ((dir == TR_FROM_UP) && (t->u1 == from))
     425                {
     426                  do_switch = TRUE;
     427                  mnew = make_new_monotone_poly(mcur, v1, v0);
     428                  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     429                  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
     430                  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     431                  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     432                }
     433              else
     434                {
     435                  mnew = make_new_monotone_poly(mcur, v0, v1);
     436                  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     437                  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     438                  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     439                  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     440                }
     441            }
     442        }
     443    }
     444  else if ((t->u0 > 0) || (t->u1 > 0)) /* no downward cusp */
     445    {
     446      if ((t->d0 > 0) && (t->d1 > 0)) /* only upward cusp */
     447        {
     448          if (_equal_to(&t->hi, &seg[t->lseg].v0))
     449            {
     450              v0 = tr[t->d1].lseg;
     451              v1 = t->lseg;
     452              retval = SP_2DN_LEFT;
     453              if (!((dir == TR_FROM_DN) && (t->d0 == from)))
     454                {
     455                  do_switch = TRUE;
     456                  mnew = make_new_monotone_poly(mcur, v1, v0);
     457                  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     458                  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     459                  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     460                  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     461                }
     462              else
     463                {
     464                  mnew = make_new_monotone_poly(mcur, v0, v1);
     465                  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     466                  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     467                  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     468                  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);           
     469                }
     470            }
     471          else
     472            {
     473              v0 = tr[t->d1].lseg;
     474              v1 = seg[t->rseg].next;
     475
     476              retval = SP_2DN_RIGHT;       
     477              if ((dir == TR_FROM_DN) && (t->d1 == from))
     478                {
     479                  do_switch = TRUE;
     480                  mnew = make_new_monotone_poly(mcur, v1, v0);
     481                  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     482                  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     483                  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     484                  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     485                }
     486              else
     487                {
     488                  mnew = make_new_monotone_poly(mcur, v0, v1);
     489                  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     490                  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     491                  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     492                  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
     493                }
     494            }
     495        }
     496      else                      /* no cusp */
     497        {
     498          if (_equal_to(&t->hi, &seg[t->lseg].v0) &&
     499              _equal_to(&t->lo, &seg[t->rseg].v0))
     500            {
     501              v0 = t->rseg;
     502              v1 = t->lseg;
     503              retval = SP_SIMPLE_LRDN;
     504              if (dir == TR_FROM_UP)
     505                {
     506                  do_switch = TRUE;
     507                  mnew = make_new_monotone_poly(mcur, v1, v0);
     508                  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     509                  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     510                  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
     511                  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     512                }
     513              else
     514                {
     515                  mnew = make_new_monotone_poly(mcur, v0, v1);
     516                  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     517                  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     518                  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     519                  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     520                }
     521            }
     522          else if (_equal_to(&t->hi, &seg[t->rseg].v1) &&
     523                   _equal_to(&t->lo, &seg[t->lseg].v1))
     524            {
     525              v0 = seg[t->rseg].next;
     526              v1 = seg[t->lseg].next;
     527
     528              retval = SP_SIMPLE_LRUP;
     529              if (dir == TR_FROM_UP)
     530                {
     531                  do_switch = TRUE;
     532                  mnew = make_new_monotone_poly(mcur, v1, v0);
     533                  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     534                  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     535                  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
     536                  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
     537                }
     538              else
     539                {
     540                  mnew = make_new_monotone_poly(mcur, v0, v1);
     541                  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
     542                  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     543                  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
     544                  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
     545                }
     546            }
     547          else                  /* no split possible */
     548            {
     549              retval = SP_NOSPLIT;
     550              traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
     551              traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
     552              traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
     553              traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);               
     554            }
     555        }
     556    }
     557
     558  return retval;
     559}
     560
     561
     562/* For each monotone polygon, find the ymax and ymin (to determine the */
     563/* two y-monotone chains) and pass on this monotone polygon for greedy */
     564/* triangulation. */
     565/* Take care not to triangulate duplicate monotone polygons */
     566
     567int triangulate_monotone_polygons(nvert, nmonpoly, op)
     568     int nvert;
     569     int nmonpoly;
     570     int op[][3];
     571{
     572  register int i;
     573  point_t ymax, ymin;
     574  int p, vfirst, posmax, posmin, v;
     575  int vcount, processed;
     576
     577#ifdef DEBUG
     578  for (i = 0; i < nmonpoly; i++)
     579    {
     580      fprintf(stderr, "\n\nPolygon %d: ", i);
     581      vfirst = mchain[mon[i]].vnum;
     582      p = mchain[mon[i]].next;
     583      fprintf (stderr, "%d ", mchain[mon[i]].vnum);
     584      while (mchain[p].vnum != vfirst)
     585        {
     586          fprintf(stderr, "%d ", mchain[p].vnum);
     587          p = mchain[p].next;
     588        }
     589    }
     590  fprintf(stderr, "\n");
     591#endif
     592
     593  op_idx = 0;
     594  for (i = 0; i < nmonpoly; i++)
     595    {
     596      vcount = 1;
     597      processed = FALSE;
     598      vfirst = mchain[mon[i]].vnum;
     599      ymax = ymin = vert[vfirst].pt;
     600      posmax = posmin = mon[i];
     601      mchain[mon[i]].marked = TRUE;
     602      p = mchain[mon[i]].next;
     603      while ((v = mchain[p].vnum) != vfirst)
     604        {
     605         if (mchain[p].marked)
     606           {
     607             processed = TRUE;
     608             break;             /* break from while */
     609           }
     610         else
     611           mchain[p].marked = TRUE;
     612
     613          if (_greater_than(&vert[v].pt, &ymax))
     614            {
     615              ymax = vert[v].pt;
     616              posmax = p;
     617            }
     618          if (_less_than(&vert[v].pt, &ymin))
     619            {
     620              ymin = vert[v].pt;
     621              posmin = p;
     622            }
     623          p = mchain[p].next;
     624          vcount++;
     625       }
     626
     627      if (processed)            /* Go to next polygon */
     628        continue;
     629     
     630      if (vcount == 3)          /* already a triangle */
     631        {
     632          op[op_idx][0] = mchain[p].vnum;
     633          op[op_idx][1] = mchain[mchain[p].next].vnum;
     634          op[op_idx][2] = mchain[mchain[p].prev].vnum;
     635          op_idx++;
     636        }
     637      else                      /* triangulate the polygon */
     638        {
     639          v = mchain[mchain[posmax].next].vnum;
     640          if (_equal_to(&vert[v].pt, &ymin))
     641            {                   /* LHS is a single line */
     642              triangulate_single_polygon(nvert, posmax, TRI_LHS, op);
     643            }
     644          else
     645            triangulate_single_polygon(nvert, posmax, TRI_RHS, op);
     646        }
     647    }
     648 
     649#ifdef DEBUG
     650  for (i = 0; i < op_idx; i++)
     651    fprintf(stderr, "tri #%d: (%d, %d, %d)\n", i, op[i][0], op[i][1],
     652           op[i][2]);
     653#endif
     654  return op_idx;
     655}
     656
     657
     658/* A greedy corner-cutting algorithm to triangulate a y-monotone
     659 * polygon in O(n) time.
     660 * Joseph O-Rourke, Computational Geometry in C.
     661 */
     662static int triangulate_single_polygon(nvert, posmax, side, op)
     663     int nvert;
     664     int posmax;
     665     int side;
     666     int op[][3];
     667{
     668  register int v;
     669  int rc[SEGSIZE], ri = 0;      /* reflex chain */
     670  int endv, tmp, vpos;
     671 
     672  if (side == TRI_RHS)          /* RHS segment is a single segment */
     673    {
     674      rc[0] = mchain[posmax].vnum;
     675      tmp = mchain[posmax].next;
     676      rc[1] = mchain[tmp].vnum;
     677      ri = 1;
     678     
     679      vpos = mchain[tmp].next;
     680      v = mchain[vpos].vnum;
     681     
     682      if ((endv = mchain[mchain[posmax].prev].vnum) == 0)
     683        endv = nvert;
     684    }
     685  else                          /* LHS is a single segment */
     686    {
     687      tmp = mchain[posmax].next;
     688      rc[0] = mchain[tmp].vnum;
     689      tmp = mchain[tmp].next;
     690      rc[1] = mchain[tmp].vnum;
     691      ri = 1;
     692
     693      vpos = mchain[tmp].next;
     694      v = mchain[vpos].vnum;
     695
     696      endv = mchain[posmax].vnum;
     697    }
     698 
     699  while ((v != endv) || (ri > 1))
     700    {
     701      if (ri > 0)               /* reflex chain is non-empty */
     702        {
     703          if (CROSS(vert[v].pt, vert[rc[ri - 1]].pt,
     704                    vert[rc[ri]].pt) > 0)
     705            {                   /* convex corner: cut if off */
     706              op[op_idx][0] = rc[ri - 1];
     707              op[op_idx][1] = rc[ri];
     708              op[op_idx][2] = v;
     709              op_idx++;     
     710              ri--;
     711            }
     712          else          /* non-convex */
     713            {           /* add v to the chain */
     714              ri++;
     715              rc[ri] = v;
     716              vpos = mchain[vpos].next;
     717              v = mchain[vpos].vnum;
     718            }
     719        }
     720      else                      /* reflex-chain empty: add v to the */
     721        {                       /* reflex chain and advance it  */
     722          rc[++ri] = v;
     723          vpos = mchain[vpos].next;
     724          v = mchain[vpos].vnum;
     725        }
     726    } /* end-while */
     727 
     728  /* reached the bottom vertex. Add in the triangle formed */
     729  op[op_idx][0] = rc[ri - 1];
     730  op[op_idx][1] = rc[ri];
     731  op[op_idx][2] = v;
     732  op_idx++;         
     733  ri--;
     734 
     735  return 0;
     736}
     737
     738
  • maptool/seidel/misc_seidel.c

     
     1#include <sys/time.h>
     2#include <math.h>
     3#include "triangulate.h"
     4
     5#ifdef __STDC__
     6extern double log2(double);
     7#else
     8extern double log2();
     9#endif
     10
     11static int choose_idx;
     12static int permute[SEGSIZE];
     13
     14
     15/* Generate a random permutation of the segments 1..n */
     16int generate_random_ordering(n)
     17     int n;
     18{
     19  struct timeval tval;
     20  struct timezone tzone;
     21  register int i;
     22  int m, st[SEGSIZE], *p;
     23 
     24  choose_idx = 1;
     25  gettimeofday(&tval, &tzone);
     26  srand48(tval.tv_sec);
     27
     28  for (i = 0; i <= n; i++)
     29    st[i] = i;
     30
     31  p = st;
     32  for (i = 1; i <= n; i++, p++)
     33    {
     34      m = lrand48() % (n + 1 - i) + 1;
     35      permute[i] = p[m];
     36      if (m != 1)
     37        p[m] = p[1];
     38    }
     39  return 0;
     40}
     41
     42 
     43/* Return the next segment in the generated random ordering of all the */
     44/* segments in S */
     45int choose_segment()
     46{
     47  int i;
     48
     49#ifdef DEBUG
     50  fprintf(stderr, "choose_segment: %d\n", permute[choose_idx]);
     51#endif
     52  return permute[choose_idx++];
     53}
     54
     55
     56#ifdef STANDALONE
     57
     58/* Read in the list of vertices from infile */
     59int read_segments(filename, genus)
     60     char *filename;
     61     int *genus;
     62{
     63  FILE *infile;
     64  int ccount;
     65  register int i;
     66  int ncontours, npoints, first, last;
     67
     68  if ((infile = fopen(filename, "r")) == NULL)
     69    {
     70      perror(filename);
     71      return -1;
     72    }
     73
     74  fscanf(infile, "%d", &ncontours);
     75  if (ncontours <= 0)
     76    return -1;
     77 
     78  /* For every contour, read in all the points for the contour. The */
     79  /* outer-most contour is read in first (points specified in */
     80  /* anti-clockwise order). Next, the inner contours are input in */
     81  /* clockwise order */
     82
     83  ccount = 0;
     84  i = 1;
     85 
     86  while (ccount < ncontours)
     87    {
     88      int j;
     89
     90      fscanf(infile, "%d", &npoints);
     91      first = i;
     92      last = first + npoints - 1;
     93      for (j = 0; j < npoints; j++, i++)
     94        {
     95          fscanf(infile, "%lf%lf", &seg[i].v0.x, &seg[i].v0.y);
     96          if (i == last)
     97            {
     98              seg[i].next = first;
     99              seg[i].prev = i-1;
     100              seg[i-1].v1 = seg[i].v0;
     101            }
     102          else if (i == first)
     103            {
     104              seg[i].next = i+1;
     105              seg[i].prev = last;
     106              seg[last].v1 = seg[i].v0;
     107            }
     108          else
     109            {
     110              seg[i].prev = i-1;
     111              seg[i].next = i+1;
     112              seg[i-1].v1 = seg[i].v0;
     113            }
     114         
     115          seg[i].is_inserted = FALSE;
     116        }
     117
     118      ccount++;
     119    }
     120
     121  *genus = ncontours - 1;
     122  return i-1;
     123}
     124
     125#endif
     126
     127
     128/* Get log*n for given n */
     129int math_logstar_n(n)
     130     int n;
     131{
     132  register int i;
     133  double v;
     134 
     135  for (i = 0, v = (double) n; v >= 1; i++)
     136    v = log2(v);
     137 
     138  return (i - 1);
     139}
     140 
     141
     142int math_N(n, h)
     143     int n;
     144     int h;
     145{
     146  register int i;
     147  double v;
     148
     149  for (i = 0, v = (int) n; i < h; i++)
     150    v = log2(v);
     151 
     152  return (int) ceil((double) 1.0*n/v);
     153}
  • maptool/seidel/construct.c

     
     1#include "triangulate.h"
     2#include <math.h>
     3
     4
     5node_t qs[QSIZE];               /* Query structure */
     6trap_t tr[TRSIZE];              /* Trapezoid structure */
     7segment_t seg[SEGSIZE];         /* Segment table */
     8
     9static int q_idx;
     10static int tr_idx;
     11
     12/* Return a new node to be added into the query tree */
     13static int newnode()
     14{
     15  if (q_idx < QSIZE)
     16    return q_idx++;
     17  else
     18    {
     19      fprintf(stderr, "newnode: Query-table overflow\n");
     20      return -1;
     21    }
     22}
     23
     24/* Return a free trapezoid */
     25static int newtrap()
     26{
     27  if (tr_idx < TRSIZE)
     28    {
     29      tr[tr_idx].lseg = -1;
     30      tr[tr_idx].rseg = -1;
     31      tr[tr_idx].state = ST_VALID;
     32      return tr_idx++;
     33    }
     34  else
     35    {
     36      fprintf(stderr, "newtrap: Trapezoid-table overflow\n");
     37      return -1;
     38    }
     39}
     40
     41
     42/* Return the maximum of the two points into the yval structure */
     43static int _max(yval, v0, v1)
     44     point_t *yval;
     45     point_t *v0;
     46     point_t *v1;
     47{
     48  if (v0->y > v1->y + C_EPS)
     49    *yval = *v0;
     50  else if (FP_EQUAL(v0->y, v1->y))
     51    {
     52      if (v0->x > v1->x + C_EPS)
     53        *yval = *v0;
     54      else
     55        *yval = *v1;
     56    }
     57  else
     58    *yval = *v1;
     59 
     60  return 0;
     61}
     62
     63
     64/* Return the minimum of the two points into the yval structure */
     65static int _min(yval, v0, v1)
     66     point_t *yval;
     67     point_t *v0;
     68     point_t *v1;
     69{
     70  if (v0->y < v1->y - C_EPS)
     71    *yval = *v0;
     72  else if (FP_EQUAL(v0->y, v1->y))
     73    {
     74      if (v0->x < v1->x)
     75        *yval = *v0;
     76      else
     77        *yval = *v1;
     78    }
     79  else
     80    *yval = *v1;
     81 
     82  return 0;
     83}
     84
     85
     86int _greater_than(v0, v1)
     87     point_t *v0;
     88     point_t *v1;
     89{
     90  if (v0->y > v1->y + C_EPS)
     91    return TRUE;
     92  else if (v0->y < v1->y - C_EPS)
     93    return FALSE;
     94  else
     95    return (v0->x > v1->x);
     96}
     97
     98
     99int _equal_to(v0, v1)
     100     point_t *v0;
     101     point_t *v1;
     102{
     103  return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
     104}
     105
     106int _greater_than_equal_to(v0, v1)
     107     point_t *v0;
     108     point_t *v1;
     109{
     110  if (v0->y > v1->y + C_EPS)
     111    return TRUE;
     112  else if (v0->y < v1->y - C_EPS)
     113    return FALSE;
     114  else
     115    return (v0->x >= v1->x);
     116}
     117
     118int _less_than(v0, v1)
     119     point_t *v0;
     120     point_t *v1;
     121{
     122  if (v0->y < v1->y - C_EPS)
     123    return TRUE;
     124  else if (v0->y > v1->y + C_EPS)
     125    return FALSE;
     126  else
     127    return (v0->x < v1->x);
     128}
     129
     130
     131/* Initilialise the query structure (Q) and the trapezoid table (T)
     132 * when the first segment is added to start the trapezoidation. The
     133 * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes
     134 *   
     135 *                4
     136 *   -----------------------------------
     137 *                \
     138 *      1          \        2
     139 *                  \
     140 *   -----------------------------------
     141 *                3
     142 */
     143
     144static int init_query_structure(segnum)
     145     int segnum;
     146{
     147  int i1, i2, i3, i4, i5, i6, i7, root;
     148  int t1, t2, t3, t4;
     149  segment_t *s = &seg[segnum];
     150
     151  q_idx = tr_idx = 1;
     152  memset((void *)tr, 0, sizeof(tr));
     153  memset((void *)qs, 0, sizeof(qs));
     154
     155  i1 = newnode();
     156  qs[i1].nodetype = T_Y;
     157  _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
     158  root = i1;
     159
     160  qs[i1].right = i2 = newnode();
     161  qs[i2].nodetype = T_SINK;
     162  qs[i2].parent = i1;
     163
     164  qs[i1].left = i3 = newnode();
     165  qs[i3].nodetype = T_Y;
     166  _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
     167  qs[i3].parent = i1;
     168 
     169  qs[i3].left = i4 = newnode();
     170  qs[i4].nodetype = T_SINK;
     171  qs[i4].parent = i3;
     172 
     173  qs[i3].right = i5 = newnode();
     174  qs[i5].nodetype = T_X;
     175  qs[i5].segnum = segnum;
     176  qs[i5].parent = i3;
     177 
     178  qs[i5].left = i6 = newnode();
     179  qs[i6].nodetype = T_SINK;
     180  qs[i6].parent = i5;
     181
     182  qs[i5].right = i7 = newnode();
     183  qs[i7].nodetype = T_SINK;
     184  qs[i7].parent = i5;
     185
     186  t1 = newtrap();               /* middle left */
     187  t2 = newtrap();               /* middle right */
     188  t3 = newtrap();               /* bottom-most */
     189  t4 = newtrap();               /* topmost */
     190
     191  tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval;
     192  tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval;
     193  tr[t4].hi.y = (double) (INFINITY);
     194  tr[t4].hi.x = (double) (INFINITY);
     195  tr[t3].lo.y = (double) -1* (INFINITY);
     196  tr[t3].lo.x = (double) -1* (INFINITY);
     197  tr[t1].rseg = tr[t2].lseg = segnum;
     198  tr[t1].u0 = tr[t2].u0 = t4;
     199  tr[t1].d0 = tr[t2].d0 = t3;
     200  tr[t4].d0 = tr[t3].u0 = t1;
     201  tr[t4].d1 = tr[t3].u1 = t2;
     202 
     203  tr[t1].sink = i6;
     204  tr[t2].sink = i7;
     205  tr[t3].sink = i4;
     206  tr[t4].sink = i2;
     207
     208  tr[t1].state = tr[t2].state = ST_VALID;
     209  tr[t3].state = tr[t4].state = ST_VALID;
     210
     211  qs[i2].trnum = t4;
     212  qs[i4].trnum = t3;
     213  qs[i6].trnum = t1;
     214  qs[i7].trnum = t2;
     215
     216  s->is_inserted = TRUE;
     217  return root;
     218}
     219
     220
     221/* Retun TRUE if the vertex v is to the left of line segment no.
     222 * segnum. Takes care of the degenerate cases when both the vertices
     223 * have the same y--cood, etc.
     224 */
     225
     226static int is_left_of(segnum, v)
     227     int segnum;
     228     point_t *v;
     229{
     230  segment_t *s = &seg[segnum];
     231  double area;
     232 
     233  if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */
     234    {
     235      if (FP_EQUAL(s->v1.y, v->y))
     236        {
     237          if (v->x < s->v1.x)
     238            area = 1.0;
     239          else
     240            area = -1.0;
     241        }
     242      else if (FP_EQUAL(s->v0.y, v->y))
     243        {
     244          if (v->x < s->v0.x)
     245            area = 1.0;
     246          else
     247            area = -1.0;
     248        }
     249      else
     250        area = CROSS(s->v0, s->v1, (*v));
     251    }
     252  else                          /* v0 > v1 */
     253    {
     254      if (FP_EQUAL(s->v1.y, v->y))
     255        {
     256          if (v->x < s->v1.x)
     257            area = 1.0;
     258          else
     259            area = -1.0;
     260        }
     261      else if (FP_EQUAL(s->v0.y, v->y))
     262        {
     263          if (v->x < s->v0.x)
     264            area = 1.0;
     265          else
     266            area = -1.0;
     267        }
     268      else
     269        area = CROSS(s->v1, s->v0, (*v));
     270    }
     271 
     272  if (area > 0.0)
     273    return TRUE;
     274  else
     275    return FALSE;
     276}
     277
     278
     279
     280/* Returns true if the corresponding endpoint of the given segment is */
     281/* already inserted into the segment tree. Use the simple test of */
     282/* whether the segment which shares this endpoint is already inserted */
     283
     284static int inserted(segnum, whichpt)
     285     int segnum;
     286     int whichpt;
     287{
     288  if (whichpt == FIRSTPT)
     289    return seg[seg[segnum].prev].is_inserted;
     290  else
     291    return seg[seg[segnum].next].is_inserted;
     292}
     293
     294/* This is query routine which determines which trapezoid does the
     295 * point v lie in. The return value is the trapezoid number.
     296 */
     297
     298int locate_endpoint(v, vo, r)
     299     point_t *v;
     300     point_t *vo;
     301     int r;
     302{
     303  node_t *rptr = &qs[r];
     304 
     305  switch (rptr->nodetype)
     306    {
     307    case T_SINK:
     308      return rptr->trnum;
     309     
     310    case T_Y:
     311      if (_greater_than(v, &rptr->yval)) /* above */
     312        return locate_endpoint(v, vo, rptr->right);
     313      else if (_equal_to(v, &rptr->yval)) /* the point is already */
     314        {                                 /* inserted. */
     315          if (_greater_than(vo, &rptr->yval)) /* above */
     316            return locate_endpoint(v, vo, rptr->right);
     317          else
     318            return locate_endpoint(v, vo, rptr->left); /* below */         
     319        }
     320      else
     321        return locate_endpoint(v, vo, rptr->left); /* below */
     322
     323    case T_X:
     324      if (_equal_to(v, &seg[rptr->segnum].v0) ||
     325               _equal_to(v, &seg[rptr->segnum].v1))
     326        {
     327          if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */
     328            {
     329              if (vo->x < v->x)
     330                return locate_endpoint(v, vo, rptr->left); /* left */
     331              else
     332                return locate_endpoint(v, vo, rptr->right); /* right */
     333            }
     334
     335          else if (is_left_of(rptr->segnum, vo))
     336            return locate_endpoint(v, vo, rptr->left); /* left */
     337          else
     338            return locate_endpoint(v, vo, rptr->right); /* right */
     339        }
     340      else if (is_left_of(rptr->segnum, v))
     341        return locate_endpoint(v, vo, rptr->left); /* left */
     342      else
     343        return locate_endpoint(v, vo, rptr->right); /* right */
     344
     345    default:
     346      fprintf(stderr, "Haggu !!!!!\n");
     347      break;
     348    }
     349}
     350
     351
     352/* Thread in the segment into the existing trapezoidation. The
     353 * limiting trapezoids are given by tfirst and tlast (which are the
     354 * trapezoids containing the two endpoints of the segment. Merges all
     355 * possible trapezoids which flank this segment and have been recently
     356 * divided because of its insertion
     357 */
     358
     359static int merge_trapezoids(segnum, tfirst, tlast, side)
     360     int segnum;
     361     int tfirst;
     362     int tlast;
     363     int side;
     364{
     365  int t, tnext, cond;
     366  int ptnext;
     367
     368  /* First merge polys on the LHS */
     369  t = tfirst;
     370  while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
     371    {
     372      if (side == S_LEFT)
     373        cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
     374                (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
     375      else
     376        cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
     377                (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));
     378     
     379      if (cond)
     380        {
     381          if ((tr[t].lseg == tr[tnext].lseg) &&
     382              (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
     383            {                                 /* merge them */
     384              /* Use the upper node as the new node i.e. t */
     385             
     386              ptnext = qs[tr[tnext].sink].parent;
     387             
     388              if (qs[ptnext].left == tr[tnext].sink)
     389                qs[ptnext].left = tr[t].sink;
     390              else
     391                qs[ptnext].right = tr[t].sink;  /* redirect parent */
     392             
     393             
     394              /* Change the upper neighbours of the lower trapezoids */
     395             
     396              if ((tr[t].d0 = tr[tnext].d0) > 0)
     397                if (tr[tr[t].d0].u0 == tnext)
     398                  tr[tr[t].d0].u0 = t;
     399                else if (tr[tr[t].d0].u1 == tnext)
     400                  tr[tr[t].d0].u1 = t;
     401             
     402              if ((tr[t].d1 = tr[tnext].d1) > 0)
     403                if (tr[tr[t].d1].u0 == tnext)
     404                  tr[tr[t].d1].u0 = t;
     405                else if (tr[tr[t].d1].u1 == tnext)
     406                  tr[tr[t].d1].u1 = t;
     407             
     408              tr[t].lo = tr[tnext].lo;
     409              tr[tnext].state = ST_INVALID; /* invalidate the lower */
     410                                            /* trapezium */
     411            }
     412          else              /* not good neighbours */
     413            t = tnext;
     414        }
     415      else                  /* do not satisfy the outer if */
     416        t = tnext;
     417     
     418    } /* end-while */
     419       
     420  return 0;
     421}
     422
     423
     424/* Add in the new segment into the trapezoidation and update Q and T
     425 * structures. First locate the two endpoints of the segment in the
     426 * Q-structure. Then start from the topmost trapezoid and go down to
     427 * the  lower trapezoid dividing all the trapezoids in between .
     428 */
     429
     430static int add_segment(segnum)
     431     int segnum;
     432{
     433  segment_t s;
     434  segment_t *so = &seg[segnum];
     435  int tu, tl, sk, tfirst, tlast, tnext;
     436  int tfirstr, tlastr, tfirstl, tlastl;
     437  int i1, i2, t, t1, t2, tn;
     438  point_t tpt;
     439  int tritop = 0, tribot = 0, is_swapped = 0;
     440  int tmptriseg;
     441
     442  s = seg[segnum];
     443  if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */
     444    {
     445      int tmp;
     446      tpt = s.v0;
     447      s.v0 = s.v1;
     448      s.v1 = tpt;
     449      tmp = s.root0;
     450      s.root0 = s.root1;
     451      s.root1 = tmp;
     452      is_swapped = TRUE;
     453    }
     454
     455  if ((is_swapped) ? !inserted(segnum, LASTPT) :
     456       !inserted(segnum, FIRSTPT))     /* insert v0 in the tree */
     457    {
     458      int tmp_d;
     459
     460      tu = locate_endpoint(&s.v0, &s.v1, s.root0);
     461      tl = newtrap();           /* tl is the new lower trapezoid */
     462      tr[tl].state = ST_VALID;
     463      tr[tl] = tr[tu];
     464      tr[tu].lo.y = tr[tl].hi.y = s.v0.y;
     465      tr[tu].lo.x = tr[tl].hi.x = s.v0.x;
     466      tr[tu].d0 = tl;     
     467      tr[tu].d1 = 0;
     468      tr[tl].u0 = tu;
     469      tr[tl].u1 = 0;
     470
     471      if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
     472        tr[tmp_d].u0 = tl;
     473      if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
     474        tr[tmp_d].u1 = tl;
     475
     476      if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
     477        tr[tmp_d].u0 = tl;
     478      if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
     479        tr[tmp_d].u1 = tl;
     480
     481      /* Now update the query structure and obtain the sinks for the */
     482      /* two trapezoids */
     483     
     484      i1 = newnode();           /* Upper trapezoid sink */
     485      i2 = newnode();           /* Lower trapezoid sink */
     486      sk = tr[tu].sink;
     487     
     488      qs[sk].nodetype = T_Y;
     489      qs[sk].yval = s.v0;
     490      qs[sk].segnum = segnum;   /* not really reqd ... maybe later */
     491      qs[sk].left = i2;
     492      qs[sk].right = i1;
     493
     494      qs[i1].nodetype = T_SINK;
     495      qs[i1].trnum = tu;
     496      qs[i1].parent = sk;
     497
     498      qs[i2].nodetype = T_SINK;
     499      qs[i2].trnum = tl;
     500      qs[i2].parent = sk;
     501
     502      tr[tu].sink = i1;
     503      tr[tl].sink = i2;
     504      tfirst = tl;
     505    }
     506  else                          /* v0 already present */
     507    {       /* Get the topmost intersecting trapezoid */
     508      tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
     509      tritop = 1;
     510    }
     511
     512
     513  if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
     514       !inserted(segnum, LASTPT))     /* insert v1 in the tree */
     515    {
     516      int tmp_d;
     517
     518      tu = locate_endpoint(&s.v1, &s.v0, s.root1);
     519
     520      tl = newtrap();           /* tl is the new lower trapezoid */
     521      tr[tl].state = ST_VALID;
     522      tr[tl] = tr[tu];
     523      tr[tu].lo.y = tr[tl].hi.y = s.v1.y;
     524      tr[tu].lo.x = tr[tl].hi.x = s.v1.x;
     525      tr[tu].d0 = tl;     
     526      tr[tu].d1 = 0;
     527      tr[tl].u0 = tu;
     528      tr[tl].u1 = 0;
     529
     530      if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
     531        tr[tmp_d].u0 = tl;
     532      if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
     533        tr[tmp_d].u1 = tl;
     534
     535      if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
     536        tr[tmp_d].u0 = tl;
     537      if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
     538        tr[tmp_d].u1 = tl;
     539     
     540      /* Now update the query structure and obtain the sinks for the */
     541      /* two trapezoids */
     542     
     543      i1 = newnode();           /* Upper trapezoid sink */
     544      i2 = newnode();           /* Lower trapezoid sink */
     545      sk = tr[tu].sink;
     546     
     547      qs[sk].nodetype = T_Y;
     548      qs[sk].yval = s.v1;
     549      qs[sk].segnum = segnum;   /* not really reqd ... maybe later */
     550      qs[sk].left = i2;
     551      qs[sk].right = i1;
     552
     553      qs[i1].nodetype = T_SINK;
     554      qs[i1].trnum = tu;
     555      qs[i1].parent = sk;
     556
     557      qs[i2].nodetype = T_SINK;
     558      qs[i2].trnum = tl;
     559      qs[i2].parent = sk;
     560
     561      tr[tu].sink = i1;
     562      tr[tl].sink = i2;
     563      tlast = tu;
     564    }
     565  else                          /* v1 already present */
     566    {       /* Get the lowermost intersecting trapezoid */
     567      tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
     568      tribot = 1;
     569    }
     570 
     571  /* Thread the segment into the query tree creating a new X-node */
     572  /* First, split all the trapezoids which are intersected by s into */
     573  /* two */
     574
     575  t = tfirst;                   /* topmost trapezoid */
     576 
     577  while ((t > 0) &&
     578         _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
     579                                /* traverse from top to bot */
     580    {
     581      int t_sav, tn_sav;
     582      sk = tr[t].sink;
     583      i1 = newnode();           /* left trapezoid sink */
     584      i2 = newnode();           /* right trapezoid sink */
     585     
     586      qs[sk].nodetype = T_X;
     587      qs[sk].segnum = segnum;
     588      qs[sk].left = i1;
     589      qs[sk].right = i2;
     590
     591      qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */
     592      qs[i1].trnum = t;
     593      qs[i1].parent = sk;
     594
     595      qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */
     596      qs[i2].trnum = tn = newtrap();
     597      tr[tn].state = ST_VALID;
     598      qs[i2].parent = sk;
     599
     600      if (t == tfirst)
     601        tfirstr = tn;
     602      if (_equal_to(&tr[t].lo, &tr[tlast].lo))
     603        tlastr = tn;
     604
     605      tr[tn] = tr[t];
     606      tr[t].sink = i1;
     607      tr[tn].sink = i2;
     608      t_sav = t;
     609      tn_sav = tn;
     610
     611      /* error */
     612
     613      if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */
     614        {
     615          fprintf(stderr, "add_segment: error\n");
     616      return 1;
     617          //break;
     618        }
     619     
     620      /* only one trapezoid below. partition t into two and make the */
     621      /* two resulting trapezoids t and tn as the upper neighbours of */
     622      /* the sole lower trapezoid */
     623     
     624      else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0))
     625        {                       /* Only one trapezoid below */
     626          if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
     627            {                   /* continuation of a chain from abv. */
     628              if (tr[t].usave > 0) /* three upper neighbours */
     629                {
     630                  if (tr[t].uside == S_LEFT)
     631                    {
     632                      tr[tn].u0 = tr[t].u1;
     633                      tr[t].u1 = -1;
     634                      tr[tn].u1 = tr[t].usave;
     635                     
     636                      tr[tr[t].u0].d0 = t;
     637                      tr[tr[tn].u0].d0 = tn;
     638                      tr[tr[tn].u1].d0 = tn;
     639                    }
     640                  else          /* intersects in the right */
     641                    {
     642                      tr[tn].u1 = -1;
     643                      tr[tn].u0 = tr[t].u1;
     644                      tr[t].u1 = tr[t].u0;
     645                      tr[t].u0 = tr[t].usave;
     646
     647                      tr[tr[t].u0].d0 = t;
     648                      tr[tr[t].u1].d0 = t;
     649                      tr[tr[tn].u0].d0 = tn;                 
     650                    }
     651                 
     652                  tr[t].usave = tr[tn].usave = 0;
     653                }
     654              else              /* No usave.... simple case */
     655                {
     656                  tr[tn].u0 = tr[t].u1;
     657                  tr[t].u1 = tr[tn].u1 = -1;
     658                  tr[tr[tn].u0].d0 = tn;
     659                }
     660            }
     661          else
     662            {                   /* fresh seg. or upward cusp */
     663              int tmp_u = tr[t].u0;
     664              int td0, td1;
     665              if (((td0 = tr[tmp_u].d0) > 0) &&
     666                  ((td1 = tr[tmp_u].d1) > 0))
     667                {               /* upward cusp */
     668                  if ((tr[td0].rseg > 0) &&
     669                      !is_left_of(tr[td0].rseg, &s.v1))
     670                    {
     671                      tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
     672                      tr[tr[tn].u0].d1 = tn;
     673                    }
     674                  else          /* cusp going leftwards */
     675                    {
     676                      tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
     677                      tr[tr[t].u0].d0 = t;
     678                    }
     679                }
     680              else              /* fresh segment */
     681                {
     682                  tr[tr[t].u0].d0 = t;
     683                  tr[tr[t].u0].d1 = tn;
     684                }             
     685            }
     686         
     687          if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
     688              FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
     689            {           /* bottom forms a triangle */
     690
     691              if (is_swapped)   
     692                tmptriseg = seg[segnum].prev;
     693              else
     694                tmptriseg = seg[segnum].next;
     695             
     696              if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
     697                {
     698                                /* L-R downward cusp */
     699                  tr[tr[t].d0].u0 = t;
     700                  tr[tn].d0 = tr[tn].d1 = -1;
     701                }
     702              else
     703                {
     704                                /* R-L downward cusp */
     705                  tr[tr[tn].d0].u1 = tn;
     706                  tr[t].d0 = tr[t].d1 = -1;
     707                }
     708            }
     709          else
     710            {
     711              if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0))
     712                {
     713                  if (tr[tr[t].d0].u0 == t) /* passes thru LHS */
     714                    {
     715                      tr[tr[t].d0].usave = tr[tr[t].d0].u1;
     716                      tr[tr[t].d0].uside = S_LEFT;
     717                    }
     718                  else
     719                    {
     720                      tr[tr[t].d0].usave = tr[tr[t].d0].u0;
     721                      tr[tr[t].d0].uside = S_RIGHT;
     722                    }               
     723                }
     724              tr[tr[t].d0].u0 = t;
     725              tr[tr[t].d0].u1 = tn;
     726            }
     727         
     728          t = tr[t].d0;
     729        }
     730
     731
     732      else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0))
     733        {                       /* Only one trapezoid below */
     734          if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
     735            {                   /* continuation of a chain from abv. */
     736              if (tr[t].usave > 0) /* three upper neighbours */
     737                {
     738                  if (tr[t].uside == S_LEFT)
     739                    {
     740                      tr[tn].u0 = tr[t].u1;
     741                      tr[t].u1 = -1;
     742                      tr[tn].u1 = tr[t].usave;
     743                     
     744                      tr[tr[t].u0].d0 = t;
     745                      tr[tr[tn].u0].d0 = tn;
     746                      tr[tr[tn].u1].d0 = tn;
     747                    }
     748                  else          /* intersects in the right */
     749                    {
     750                      tr[tn].u1 = -1;
     751                      tr[tn].u0 = tr[t].u1;
     752                      tr[t].u1 = tr[t].u0;
     753                      tr[t].u0 = tr[t].usave;
     754
     755                      tr[tr[t].u0].d0 = t;
     756                      tr[tr[t].u1].d0 = t;
     757                      tr[tr[tn].u0].d0 = tn;                 
     758                    }
     759                 
     760                  tr[t].usave = tr[tn].usave = 0;
     761                }
     762              else              /* No usave.... simple case */
     763                {
     764                  tr[tn].u0 = tr[t].u1;
     765                  tr[t].u1 = tr[tn].u1 = -1;
     766                  tr[tr[tn].u0].d0 = tn;
     767                }
     768            }
     769          else
     770            {                   /* fresh seg. or upward cusp */
     771              int tmp_u = tr[t].u0;
     772              int td0, td1;
     773              if (((td0 = tr[tmp_u].d0) > 0) &&
     774                  ((td1 = tr[tmp_u].d1) > 0))
     775                {               /* upward cusp */
     776                  if ((tr[td0].rseg > 0) &&
     777                      !is_left_of(tr[td0].rseg, &s.v1))
     778                    {
     779                      tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
     780                      tr[tr[tn].u0].d1 = tn;
     781                    }
     782                  else
     783                    {
     784                      tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
     785                      tr[tr[t].u0].d0 = t;
     786                    }
     787                }
     788              else              /* fresh segment */
     789                {
     790                  tr[tr[t].u0].d0 = t;
     791                  tr[tr[t].u0].d1 = tn;
     792                }
     793            }
     794         
     795          if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
     796              FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
     797            {           /* bottom forms a triangle */
     798              int tmpseg;
     799
     800              if (is_swapped)   
     801                tmptriseg = seg[segnum].prev;
     802              else
     803                tmptriseg = seg[segnum].next;
     804
     805              if ((tmpseg > 0) && is_left_of(tmpseg, &s.v0))
     806                {
     807                  /* L-R downward cusp */
     808                  tr[tr[t].d1].u0 = t;
     809                  tr[tn].d0 = tr[tn].d1 = -1;
     810                }
     811              else
     812                {
     813                  /* R-L downward cusp */
     814                  tr[tr[tn].d1].u1 = tn;
     815                  tr[t].d0 = tr[t].d1 = -1;
     816                }
     817            }           
     818          else
     819            {
     820              if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0))
     821                {
     822                  if (tr[tr[t].d1].u0 == t) /* passes thru LHS */
     823                    {
     824                      tr[tr[t].d1].usave = tr[tr[t].d1].u1;
     825                      tr[tr[t].d1].uside = S_LEFT;
     826                    }
     827                  else
     828                    {
     829                      tr[tr[t].d1].usave = tr[tr[t].d1].u0;
     830                      tr[tr[t].d1].uside = S_RIGHT;
     831                    }               
     832                }
     833              tr[tr[t].d1].u0 = t;
     834              tr[tr[t].d1].u1 = tn;
     835            }
     836         
     837          t = tr[t].d1;
     838        }
     839
     840      /* two trapezoids below. Find out which one is intersected by */
     841      /* this segment and proceed down that one */
     842     
     843      else
     844        {
     845          int tmpseg = tr[tr[t].d0].rseg;
     846          double y0, yt;
     847          point_t tmppt;
     848          int tnext, i_d0, i_d1;
     849
     850          i_d0 = i_d1 = FALSE;
     851          if (FP_EQUAL(tr[t].lo.y, s.v0.y))
     852            {
     853              if (tr[t].lo.x > s.v0.x)
     854                i_d0 = TRUE;
     855              else
     856                i_d1 = TRUE;
     857            }
     858          else
     859            {
     860              tmppt.y = y0 = tr[t].lo.y;
     861              yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y);
     862              tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);
     863             
     864              if (_less_than(&tmppt, &tr[t].lo))
     865                i_d0 = TRUE;
     866              else
     867                i_d1 = TRUE;
     868            }
     869         
     870          /* check continuity from the top so that the lower-neighbour */
     871          /* values are properly filled for the upper trapezoid */
     872
     873          if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
     874            {                   /* continuation of a chain from abv. */
     875              if (tr[t].usave > 0) /* three upper neighbours */
     876                {
     877                  if (tr[t].uside == S_LEFT)
     878                    {
     879                      tr[tn].u0 = tr[t].u1;
     880                      tr[t].u1 = -1;
     881                      tr[tn].u1 = tr[t].usave;
     882                     
     883                      tr[tr[t].u0].d0 = t;
     884                      tr[tr[tn].u0].d0 = tn;
     885                      tr[tr[tn].u1].d0 = tn;
     886                    }
     887                  else          /* intersects in the right */
     888                    {
     889                      tr[tn].u1 = -1;
     890                      tr[tn].u0 = tr[t].u1;
     891                      tr[t].u1 = tr[t].u0;
     892                      tr[t].u0 = tr[t].usave;
     893
     894                      tr[tr[t].u0].d0 = t;
     895                      tr[tr[t].u1].d0 = t;
     896                      tr[tr[tn].u0].d0 = tn;                 
     897                    }
     898                 
     899                  tr[t].usave = tr[tn].usave = 0;
     900                }
     901              else              /* No usave.... simple case */
     902                {
     903                  tr[tn].u0 = tr[t].u1;
     904                  tr[tn].u1 = -1;
     905                  tr[t].u1 = -1;
     906                  tr[tr[tn].u0].d0 = tn;
     907                }
     908            }
     909          else
     910            {                   /* fresh seg. or upward cusp */
     911              int tmp_u = tr[t].u0;
     912              int td0, td1;
     913              if (((td0 = tr[tmp_u].d0) > 0) &&
     914                  ((td1 = tr[tmp_u].d1) > 0))
     915                {               /* upward cusp */
     916                  if ((tr[td0].rseg > 0) &&
     917                      !is_left_of(tr[td0].rseg, &s.v1))
     918                    {
     919                      tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
     920                      tr[tr[tn].u0].d1 = tn;
     921                    }
     922                  else
     923                    {
     924                      tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
     925                      tr[tr[t].u0].d0 = t;
     926                    }
     927                }
     928              else              /* fresh segment */
     929                {
     930                  tr[tr[t].u0].d0 = t;
     931                  tr[tr[t].u0].d1 = tn;
     932                }
     933            }
     934         
     935          if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
     936              FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
     937            {
     938              /* this case arises only at the lowest trapezoid.. i.e.
     939                 tlast, if the lower endpoint of the segment is
     940                 already inserted in the structure */
     941             
     942              tr[tr[t].d0].u0 = t;
     943              tr[tr[t].d0].u1 = -1;
     944              tr[tr[t].d1].u0 = tn;
     945              tr[tr[t].d1].u1 = -1;
     946
     947              tr[tn].d0 = tr[t].d1;
     948              tr[t].d1 = tr[tn].d1 = -1;
     949             
     950              tnext = tr[t].d1;       
     951            }
     952          else if (i_d0)
     953                                /* intersecting d0 */
     954            {
     955              tr[tr[t].d0].u0 = t;
     956              tr[tr[t].d0].u1 = tn;
     957              tr[tr[t].d1].u0 = tn;
     958              tr[tr[t].d1].u1 = -1;
     959             
     960              /* new code to determine the bottom neighbours of the */
     961              /* newly partitioned trapezoid */
     962             
     963              tr[t].d1 = -1;
     964
     965              tnext = tr[t].d0;
     966            }
     967          else                  /* intersecting d1 */
     968            {
     969              tr[tr[t].d0].u0 = t;
     970              tr[tr[t].d0].u1 = -1;
     971              tr[tr[t].d1].u0 = t;
     972              tr[tr[t].d1].u1 = tn;
     973
     974              /* new code to determine the bottom neighbours of the */
     975              /* newly partitioned trapezoid */
     976             
     977              tr[tn].d0 = tr[t].d1;
     978              tr[tn].d1 = -1;
     979             
     980              tnext = tr[t].d1;
     981            }       
     982         
     983          t = tnext;
     984        }
     985     
     986      tr[t_sav].rseg = tr[tn_sav].lseg  = segnum;
     987    } /* end-while */
     988 
     989  /* Now combine those trapezoids which share common segments. We can */
     990  /* use the pointers to the parent to connect these together. This */
     991  /* works only because all these new trapezoids have been formed */
     992  /* due to splitting by the segment, and hence have only one parent */
     993
     994  tfirstl = tfirst;
     995  tlastl = tlast;
     996  merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
     997  merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);
     998
     999  seg[segnum].is_inserted = TRUE;
     1000  return 0;
     1001}
     1002
     1003
     1004/* Update the roots stored for each of the endpoints of the segment.
     1005 * This is done to speed up the location-query for the endpoint when
     1006 * the segment is inserted into the trapezoidation subsequently
     1007 */
     1008static int find_new_roots(segnum)
     1009     int segnum;
     1010{
     1011  segment_t *s = &seg[segnum];
     1012 
     1013  if (s->is_inserted)
     1014    return 0;
     1015
     1016  s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
     1017  s->root0 = tr[s->root0].sink;
     1018
     1019  s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
     1020  s->root1 = tr[s->root1].sink; 
     1021  return 0;
     1022}
     1023
     1024
     1025/* Main routine to perform trapezoidation */
     1026int construct_trapezoids(nseg)
     1027     int nseg;
     1028{
     1029  register int i;
     1030  int root, h;
     1031 
     1032  /* Add the first segment and get the query structure and trapezoid */
     1033  /* list initialised */
     1034
     1035  root = init_query_structure(choose_segment());
     1036
     1037  for (i = 1; i <= nseg; i++)
     1038    seg[i].root0 = seg[i].root1 = root;
     1039 
     1040  for (h = 1; h <= math_logstar_n(nseg); h++)
     1041    {
     1042      for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++)
     1043            if(add_segment(choose_segment())) {
     1044                  return 1;
     1045            }
     1046     
     1047      /* Find a new root for each of the segment endpoints */
     1048      for (i = 1; i <= nseg; i++)
     1049        find_new_roots(i);
     1050    }
     1051 
     1052  for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
     1053    if(add_segment(choose_segment())) {
     1054      return 1;
     1055    }
     1056
     1057  return 0;
     1058}
     1059
     1060
  • maptool/seidel/triangulate.h

     
     1#ifndef _triangulate_h
     2#define _triangulate_h
     3
     4#include <sys/types.h>
     5#include <stdlib.h>
     6#include <stdio.h>
     7#include "types.h"
     8
     9typedef struct {
     10  double x, y;
     11} point_t, vector_t;
     12
     13
     14/* Segment attributes */
     15
     16typedef struct {       
     17  point_t v0, v1;               /* two endpoints */
     18  int is_inserted;              /* inserted in trapezoidation yet ? */
     19  int root0, root1;             /* root nodes in Q */
     20  int next;                     /* Next logical segment */
     21  int prev;                     /* Previous segment */
     22} segment_t;
     23
     24
     25/* Trapezoid attributes */
     26
     27typedef struct {
     28  int lseg, rseg;               /* two adjoining segments */
     29  point_t hi, lo;               /* max/min y-values */
     30  int u0, u1;
     31  int d0, d1;
     32  int sink;                     /* pointer to corresponding in Q */
     33  int usave, uside;             /* I forgot what this means */
     34  int state;
     35} trap_t;
     36
     37
     38/* Node attributes for every node in the query structure */
     39
     40typedef struct {
     41  int nodetype;                 /* Y-node or S-node */
     42  int segnum;
     43  point_t yval;
     44  int trnum;
     45  int parent;                   /* doubly linked DAG */
     46  int left, right;              /* children */
     47} node_t;
     48
     49
     50typedef struct {
     51  int vnum;
     52  int next;                     /* Circularly linked list  */
     53  int prev;                     /* describing the monotone */
     54  int marked;                   /* polygon */
     55} monchain_t;                   
     56
     57
     58typedef struct {
     59  point_t pt;
     60  int vnext[4];                 /* next vertices for the 4 chains */
     61  int vpos[4];                  /* position of v in the 4 chains */
     62  int nextfree;
     63} vertexchain_t;
     64
     65
     66/* Node types */
     67
     68#define T_X     1
     69#define T_Y     2
     70#define T_SINK  3
     71
     72
     73#define SEGSIZE 1000000         /* max# of segments. Determines how */
     74                                /* many points can be specified as */
     75                                /* input. If your datasets have large */
     76                                /* number of points, increase this */
     77                                /* value accordingly. */
     78
     79#define QSIZE   8*SEGSIZE       /* maximum table sizes */
     80#define TRSIZE  4*SEGSIZE       /* max# trapezoids */
     81
     82#define TRUE  1
     83#define FALSE 0
     84
     85#define FIRSTPT 1               /* checking whether pt. is inserted */
     86#define LASTPT  2
     87
     88
     89#define INFINITY 1<<30
     90//#define C_EPS 1.0e-7          /* tolerance value: Used for making */
     91#define C_EPS 1.0e-7            /* tolerance value: Used for making */
     92                                /* all decisions about collinearity or */
     93                                /* left/right of segment. Decrease */
     94                                /* this value if the input points are */
     95                                /* spaced very close together */
     96
     97
     98#define S_LEFT 1                /* for merge-direction */
     99#define S_RIGHT 2
     100
     101
     102#define ST_VALID 1              /* for trapezium state */
     103#define ST_INVALID 2
     104
     105
     106#define SP_SIMPLE_LRUP 1        /* for splitting trapezoids */
     107#define SP_SIMPLE_LRDN 2
     108#define SP_2UP_2DN     3
     109#define SP_2UP_LEFT    4
     110#define SP_2UP_RIGHT   5
     111#define SP_2DN_LEFT    6
     112#define SP_2DN_RIGHT   7
     113#define SP_NOSPLIT    -1       
     114
     115#define TR_FROM_UP 1            /* for traverse-direction */
     116#define TR_FROM_DN 2
     117
     118#define TRI_LHS 1
     119#define TRI_RHS 2
     120
     121
     122#define MAX(a, b) (((a) > (b)) ? (a) : (b))
     123#define MIN(a, b) (((a) < (b)) ? (a) : (b))
     124
     125#define CROSS(v0, v1, v2) (((v1).x - (v0).x)*((v2).y - (v0).y) - \
     126                           ((v1).y - (v0).y)*((v2).x - (v0).x))
     127
     128#define DOT(v0, v1) ((v0).x * (v1).x + (v0).y * (v1).y)
     129
     130#define FP_EQUAL(s, t) (fabs(s - t) <= C_EPS)
     131
     132
     133
     134/* Global variables */
     135
     136extern node_t qs[QSIZE];                /* Query structure */
     137extern trap_t tr[TRSIZE];               /* Trapezoid structure */
     138extern segment_t seg[SEGSIZE];          /* Segment table */
     139
     140
     141/* Functions */
     142
     143extern int monotonate_trapezoids(int);
     144extern int triangulate_monotone_polygons(int, int, int (*)[3]);
     145
     146extern int _greater_than(point_t *, point_t *);
     147extern int _equal_to(point_t *, point_t *);
     148extern int _greater_than_equal_to(point_t *, point_t *);
     149extern int _less_than(point_t *, point_t *);
     150extern int locate_endpoint(point_t *, point_t *, int);
     151extern int construct_trapezoids(int);
     152
     153extern int generate_random_ordering(int);
     154extern int choose_segment(void);
     155extern int read_segments(char *, int *);
     156extern int math_logstar_n(int);
     157extern int math_N(int, int);
     158
     159#endif /* triangulate_h */
  • maptool/seidel/interface.h

     
     1#ifndef __interface_h
     2#define __interface_h
     3
     4#include "types.h"
     5
     6#define TRUE 1
     7#define FALSE 0
     8
     9extern int triangulate_polygon(s32, s32*, s32 (*)[2], s32 (*)[3]);
     10extern int is_point_inside_polygon(double *);
     11
     12#endif /* __interface_h */
  • maptool/geom.c

     
    1919#include <string.h>
    2020#include <math.h>
    2121#include "maptool.h"
     22#include "transform.h"
    2223
    2324void
    2425geom_coord_copy(struct coord *from, struct coord *to, int count, int reverse)
     
    357358        return ret;
    358359}
    359360
     361/*
     362 * create groups of segments by merging segments with common endpoint
     363 */
     364GList* geom_poly_segments_group(GList *in, GList *out)
     365{
     366        out = g_list_copy(in);
     367
     368        int changed;
     369        do {
     370                changed = 0;
     371                GList*l1 = out;
     372                while(l1) {
     373                        struct geom_poly_segment *l1_ = l1->data;
     374                        int size_1 = (l1_->last-l1_->first)+1;
     375                        if(((l1_->type!=geom_poly_segment_type_way_outer) && (l1_->type!=geom_poly_segment_type_way_inner)) || (size_1<1)) {
     376                                GList*l1_d = l1;
     377                                l1 = g_list_next(l1);
     378                                if(!l1) {
     379                                        break;
     380                                }
     381                                l1_ = l1->data;
     382                                out = g_list_remove(out, l1_d->data);
     383                                continue;
     384                        }
     385                        GList*l2 = l1;
     386                        while(l2) {
     387                                struct geom_poly_segment *l2_ = l2->data;
     388                                int size_2 = (l2_->last-l2_->first)+1;
     389                                if(((l2_->type!=geom_poly_segment_type_way_outer) && (l2_->type!=geom_poly_segment_type_way_inner)) || (size_2<1)) {
     390                                        GList*l2_d = l2;
     391                                        l2 = g_list_next(l2);
     392                                        if(!l2) {
     393                                                break;
     394                                        }
     395                                        l2_ = l2->data;
     396                                        out = g_list_remove(out, l2_d->data);
     397                                        continue;
     398                                }
     399                                //suppose that inner segments are only connected with other inner segments
     400                                if( l1!=l2 && l1_ && l2_ && l1_->first && l2_->last &&
     401                                        (coord_equal(l1_->first, l2_->first) || /* have_common_endpoints(l1,l2)*/
     402                                                coord_equal(l1_->first, l2_->last) ||
     403                                                coord_equal(l1_->last, l2_->first) ||
     404                                                coord_equal(l1_->last, l2_->last))              &&
     405                                        !coord_equal(l1_->first,l1_->last) && !coord_equal(l2_->first,l2_->last) /*do not group closed polys*/ ) {
     406                                        int reverse_1 = 0;
     407                                        int reverse_2 = 0;
     408                                        int order_12 = 1;
     409                                        //merge coords of (l1,l2) (append coords to the front or rear, perhaps reversed)
     410                                        //l1.front==l2.front
     411                                                //l2 reversed ... l1 inorder
     412                                        if (coord_equal(l1_->first,l2_->first)) {
     413                                                reverse_1 = 0;
     414                                                reverse_2 = 1;
     415                                                order_12  = 0;
     416                                        }
     417                                        //l1.rear==l2.front
     418                                                //l1 inorder ... l2 inorder
     419                                        else if (coord_equal(l1_->last,l2_->first)) {
     420                                                reverse_1 = 0;
     421                                                reverse_2 = 0;
     422                                                order_12  = 1;
     423                                        }
     424                                        //l1.front==l2.rear
     425                                                //l2 inorder ... l1 inorder
     426                                        else if (coord_equal(l1_->first,l2_->last)) {
     427                                                reverse_1 = 0;
     428                                                reverse_2 = 0;
     429                                                order_12  = 0;
     430                                        }
     431                                                        //l1.rear==l2.rear
     432                                                //l1 inorder ... l2 reversed
     433                                        else if (coord_equal(l1_->last,l2_->last)) {
     434                                                reverse_1 = 0;
     435                                                reverse_2 = 1;
     436                                                order_12  = 1;
     437                                        }
     438                                        //reverse l1 or l2 if necessary
     439                                        if (reverse_1) {
     440                                                int ii;
     441                                                for (ii=0;ii<size_1/2;++ii) {
     442                                                        struct coord tmp = l1_->first[ii];
     443                                                        l1_->first[ii] = l1_->first[size_1-ii-1];
     444                                                        l1_->first[size_1-ii-1] = tmp;
     445                                                }
     446                                        }
     447                                        if (reverse_2) {
     448                                                int ii;
     449                                                for (ii=0;ii<size_2/2;++ii) {
     450                                                        struct coord tmp = l2_->first[ii];
     451                                                        l2_->first[ii] = l2_->first[size_2-ii-1];
     452                                                        l2_->first[size_2-ii-1] = tmp;
     453                                                }
     454                                        }
     455                                        //append l1 to l2 or l2 to l1 appropriately (with the common point added only once)
     456                                        if(order_12 && size_2>1) {
     457                                                l1_->first = g_realloc(l1_->first,sizeof(struct coord)*(size_1+size_2)-1);
     458                                                memcpy(l1_->first+size_1,l2_->first,(size_2-1)*sizeof(struct coord));
     459                                                l1_->last += size_2-1;
     460                                        } else if( !order_12 && size_1>1) {
     461                                                l2_->first = g_realloc(l2_->first,sizeof(struct coord)*(size_1+size_2-1));
     462                                                memcpy(l2_->first+size_2,l1_->first,(size_1-1)*sizeof(struct coord));
     463                                                l2_->last += size_1-1;
     464                                        }
     465                                        //remove l1 or l2 appropriately
     466                                        if (order_12) {
     467                                                out = g_list_remove(out, l2->data);
     468                                        } else {
     469                                                out = g_list_remove(out, l1->data);
     470                                        }
     471                                        changed=1;
     472                                }
     473                                l2 = g_list_next(l2);
     474                        }
     475                        l1 = g_list_next(l1);
     476                }
     477        } while(changed);
     478return out;
     479}
     480
     481
    360482int
    361483geom_poly_segments_point_inside(GList *in, struct coord *c)
    362484{
     
    620742                item_bin_write_clipped(ib_in, param, out);
    621743        }
    622744}
     745
     746
     747
     748
     749static int coord_dot_prod (struct coord* c1, struct coord* c2)
     750{
     751        return c1->x*c2->x + c1->y*c2->y;
     752}
     753
     754
     755static void coord_sub (struct coord* c1, struct coord* c2, struct coord* cout)
     756{
     757        cout->x = c1->x - c2->x;
     758        cout->y = c1->y - c2->y;
     759}
     760
     761double DistanceSquared(struct coord* v0, struct coord* v1, struct coord*p)
     762{
     763
     764/*handle the case of vertical line segments*/
     765        if (v0->x==v1->x && p->x==v0->x &&
     766                ( (v0->y<=p->y && p->y<=v1->y) || (v1->y<=p->y && p->y<=v0->y) )
     767        ) {
     768                return 0.0;
     769        }
     770
     771        if (v0->y==v1->y && p->y==v0->y &&
     772                ( (v0->x<=p->x && p->x<=v1->x) || (v1->x<=p->x && p->x<=v0->x) )
     773        ) {
     774                return 0.0;
     775        }
     776
     777
     778        double vx = v0->x - p->x, vy = v0->y - p->y, ux = v1->x - v0->x, uy = v1->y - v0->y;
     779        double length = ux * ux + uy * uy;
     780
     781        double det = (-vx * ux) + (-vy * uy); //if this is < 0 or > length then its outside the line segment
     782        if (det < 0) {
     783                return (v0->x - p->x) * (v0->x - p->x) + (v0->y - p->y) * (v0->y - p->y);
     784        }
     785        if (det > length) {
     786                return (v1->x - p->x) * (v1->x - p->x) + (v1->y - p->y) * (v1->y - p->y);
     787        }
     788
     789        det = ux * vy - uy * vx;
     790        return (det * det) / length;
     791}
     792
     793
     794//TODO merge intersect tests
     795int 
     796self_intersect_test_segment(struct geom_poly_segment *gs)
     797{
     798        int vertices_count = (gs->last-gs->first)+1;
     799        if(vertices_count<4) {
     800                return 0;
     801        }
     802        struct coord*vertices = (struct coord*)gs->first;
     803        if(coord_equal(&vertices[0],&vertices[vertices_count-1])) {
     804                vertices_count--;
     805        }
     806        int i;
     807        for (i = 0; i < vertices_count; ++i) {
     808                if (i < vertices_count - 1) {
     809                    int h;
     810                    for (h = i + 1; h < vertices_count; ++h)
     811                    {
     812                        // Do two vertices lie on top of one another?
     813                        if ( i!=0 && h-i!=1 && h-i!=-1 && vertices[i].x == vertices[h].x && vertices[i].y == vertices[h].y ) {
     814                            return 1;
     815                        }
     816                    }
     817                }
     818
     819                int j = (i + 1) % vertices_count;
     820                struct coord iToj;
     821                coord_sub(&vertices[j], &vertices[i], &iToj);
     822                struct coord iTojNormal;
     823                iTojNormal.x =  iToj.y;
     824                iTojNormal.y = -iToj.x;
     825                // i is the first vertex and j is the second
     826                int startK = (j + 1) % vertices_count;
     827                int endK = (i - 1 + vertices_count) % vertices_count;
     828                endK += startK < endK ? 0 : startK + 1;
     829                int k = startK;
     830                struct coord iTok;
     831                coord_sub(&vertices[k], &vertices[i], &iTok);
     832                int onLeftSide = coord_dot_prod(&iTok, &iTojNormal) >= 0;
     833                struct coord prevK = vertices[k];
     834                ++k;
     835                for (; k <= endK; ++k)
     836                {
     837                    int modK = k % vertices_count;
     838                    coord_sub(&vertices[modK], &vertices[i], &iTok);
     839                    if (onLeftSide != coord_dot_prod(&iTok, &iTojNormal) >= 0)
     840                    {
     841                        struct coord prevKtoK, idiff, jdiff;
     842                        coord_sub(&vertices[modK], &prevK, &prevKtoK);
     843                        struct coord prevKtoKNormal;
     844                        prevKtoKNormal.x = prevKtoK.y;
     845                        prevKtoKNormal.y = -prevKtoK.x;
     846                        coord_sub(&vertices[i],&prevK,&idiff);
     847                        coord_sub(&vertices[j],&prevK,&jdiff);
     848                        if ((coord_dot_prod(&idiff, &prevKtoKNormal) >= 0) != (coord_dot_prod(&jdiff, &prevKtoKNormal) >= 0))
     849                        {
     850                               return 1;
     851                        }
     852                    }
     853                    onLeftSide = coord_dot_prod(&iTok, &iTojNormal) > 0;
     854                    prevK = vertices[modK];
     855                }
     856            }
     857                        //check if a vertex lies on a lineseg defined by other vertices
     858                        struct coord* c = gs->first;
     859                        struct coord* c2;
     860                        while (c+1<=gs->last) {
     861                                c2 = gs->first;
     862                                while (c2<=gs->last) {
     863                                        if ( !coord_equal(c,c2) && !coord_equal(c+1,c2) && DistanceSquared(c,c+1,c2)<.01) {
     864                                                return 1;
     865                                        }
     866                                        c2++;
     867                                }
     868                                c++;
     869                        }
     870
     871                        if (!coord_equal(gs->first,gs->last)) {
     872                                c2 = gs->first;
     873                                while (c2<=gs->last) {
     874                                        if ( !coord_equal(gs->first,c2) && !coord_equal(gs->last,c2) && DistanceSquared(gs->first,gs->last,c2)<.01) {
     875                                                        return 1;
     876                                        }
     877                                        c2++;
     878                                }
     879                        }
     880            return 0;
     881
     882}
     883
     884
     885int self_intersect_test(struct item_bin*ib)
     886{
     887        int vertices_count = ib->clen/2;
     888        if(vertices_count<4) {
     889                return 0;
     890        }
     891        struct coord*vertices = (struct coord*)(ib+1);
     892
     893        if(coord_equal(&vertices[0],&vertices[vertices_count-1])) {
     894                vertices_count--;
     895        }
     896        int i;
     897        for (i = 0; i < vertices_count; ++i) {
     898                if (i < vertices_count - 1) {
     899                    int h;
     900                    for (h = i + 1; h < vertices_count; ++h)
     901                    {
     902                        // Do two vertices lie on top of one another?
     903                        if ( i!=0 && h-i!=1 && h-i!=-1 && vertices[i].x == vertices[h].x && vertices[i].y == vertices[h].y ) {
     904                            return 1;
     905                        }
     906                    }
     907                }
     908
     909                int j = (i + 1) % vertices_count;
     910                struct coord iToj;
     911                coord_sub(&vertices[j], &vertices[i], &iToj);
     912                struct coord iTojNormal;
     913                iTojNormal.x =  iToj.y;
     914                iTojNormal.y = -iToj.x;
     915                // i is the first vertex and j is the second
     916                int startK = (j + 1) % vertices_count;
     917                int endK = (i - 1 + vertices_count) % vertices_count;
     918                endK += startK < endK ? 0 : startK + 1;
     919                int k = startK;
     920                struct coord iTok;
     921                coord_sub(&vertices[k], &vertices[i], &iTok);
     922                int onLeftSide = coord_dot_prod(&iTok, &iTojNormal) >= 0;
     923                struct coord prevK = vertices[k];
     924                ++k;
     925                for (; k <= endK; ++k)
     926                {
     927                    int modK = k % vertices_count;
     928                    coord_sub(&vertices[modK], &vertices[i], &iTok);
     929                    if (onLeftSide != coord_dot_prod(&iTok, &iTojNormal) >= 0)
     930                    {
     931                        struct coord prevKtoK, idiff, jdiff;
     932                        coord_sub(&vertices[modK], &prevK, &prevKtoK);
     933                        struct coord prevKtoKNormal;
     934                        prevKtoKNormal.x = prevKtoK.y;
     935                        prevKtoKNormal.y = -prevKtoK.x;
     936                        coord_sub(&vertices[i],&prevK,&idiff);
     937                        coord_sub(&vertices[j],&prevK,&jdiff);
     938                        if ((coord_dot_prod(&idiff, &prevKtoKNormal) >= 0) != (coord_dot_prod(&jdiff, &prevKtoKNormal) >= 0))
     939                        {
     940                            if(1/*i!=j-1 && i!=j+1*/) {
     941                               return 1;
     942                            }
     943                        }
     944                    }
     945                    onLeftSide = coord_dot_prod(&iTok, &iTojNormal) > 0;
     946                    prevK = vertices[modK];
     947                }
     948            }
     949            return 0;
     950}
     951
     952/*
     953 * returns true if coord c is inside poly segment poly
     954 */
     955int geom_point_in_polygon(struct coord*c, struct geom_poly_segment*poly)
     956{
     957        //TODO handle closed/opened polygon cases
     958        int counter = 0;
     959        int i,N;
     960        double xinters;
     961        struct coord p1,p2;
     962
     963        if(!c || !poly) {
     964                return 0;
     965        }
     966        N = poly->last-poly->first+1;
     967        p1 = poly->first[0];
     968       
     969        for (i=1;i<=N;i++) {
     970                p2 = poly->first[i % N];
     971                if (c->y > MIN(p1.y,p2.y)) {
     972                        if (c->y <= MAX(p1.y,p2.y)) {
     973                                if (c->x <= MAX(p1.x,p2.x)) {
     974                                        if (p1.y != p2.y) {
     975                                                xinters = (c->y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
     976                                                if (p1.x == p2.x || c->x <= xinters)
     977                                                        counter++;
     978                                        }
     979                                }
     980                        }
     981                }
     982                p1 = p2;
     983        }
     984
     985        if (counter % 2 == 0)
     986                return 0;
     987        else
     988                return 1;
     989}
     990
     991/*
     992 * returns 1 if poly1 is approximately inside poly2 (every point of p1 is inside poly2)
     993 */
     994int geom_polygon_in_polygon(struct geom_poly_segment*poly1, struct geom_poly_segment*poly2)
     995{
     996        int ret=1;
     997        struct coord *c = poly1->first;
     998        do {
     999                if(!geom_point_in_polygon(c,poly2)) {
     1000                        ret=0;
     1001                        break;
     1002                }
     1003        } while(c++ != poly1->last);
     1004        return ret;
     1005}
     1006
     1007
     1008int geom_polygon_is_clockwise(struct geom_poly_segment*poly)
     1009{
     1010        double sum=0;   //sum of signed areas
     1011        struct coord *c = poly->first;
     1012        struct coord *c2;
     1013        do {
     1014                c2 = c+1;
     1015                sum += (double)c->x*c2->y-(double)c2->x*c->y;
     1016        } while(++c <= (poly->last-1));
     1017        if(!coord_equal(poly->first,poly->last)) {
     1018                c =poly->last;
     1019                c2=poly->first;
     1020                sum += (double)c->x*c2->y-(double)c2->x*c->y;
     1021        }
     1022        return sum<0 ? 1 : 0;
     1023}
     1024
     1025void geom_polygon_reverse(struct geom_poly_segment*poly)
     1026{
     1027        int n = poly->last-poly->first+1;
     1028        int ii;
     1029        for(ii=0;ii<n/2;++ii) {
     1030                struct coord tmp = poly->first[ii];
     1031                poly->first[ii] = poly->first[n-1-ii];
     1032                poly->first[n-1-ii] = tmp;
     1033        }
     1034}
     1035
     1036/*
     1037 * removes subsequent equivalent coords from segment
     1038 */
     1039void geom_polygon_unique(struct geom_poly_segment*poly)
     1040{
     1041        struct coord *c  = poly->first;
     1042        struct coord *dst = NULL;
     1043        int unique_cnt=2;
     1044        int rearrange=0;
     1045        if(poly->last<=poly->first) {
     1046                return;
     1047        }
     1048        do {
     1049                struct coord *c2 = c+1;
     1050                if(!coord_equal(c,c2)) {
     1051                        ++unique_cnt;
     1052                } else {
     1053                        rearrange=1;
     1054                }
     1055        } while(++c != poly->last);
     1056
     1057        if(rearrange) {
     1058                c = poly->first;
     1059                struct coord*new_segment_data = g_malloc(sizeof(struct coord)*unique_cnt);
     1060                dst=new_segment_data;
     1061                do {
     1062                        struct coord *c2 = c+1;
     1063                        *dst++ = *c;
     1064                        while(c2<=poly->last && coord_equal(c,c2)) {
     1065                                ++c;
     1066                                ++c2;
     1067                        }
     1068                } while(++c <= poly->last);
     1069                g_free(poly->first);
     1070                poly->first=new_segment_data;
     1071                poly->last=dst-1;
     1072        }
     1073}
     1074
     1075int geom_triangle_ccw(struct coord*a,struct coord*b,struct coord*c)
     1076{
     1077        return (c->y-a->y)*(b->x-a->x) >= (b->y-a->y)*(c->x-a->x);
     1078}
     1079
     1080int geom_lineseg_intersect(struct coord*a,struct coord*b,struct coord*c,struct coord*d)
     1081{
     1082        if(DistanceSquared(a,b,c)<.1 || DistanceSquared(a,b,d)<.1 ||
     1083        DistanceSquared(c,d,a)<.1 || DistanceSquared(c,d,b)<.1  ) {
     1084                return 1;
     1085        }
     1086
     1087        int d1 = transform_distance_line_sq(a,b,c,NULL);
     1088        int d2 = transform_distance_line_sq(a,b,d,NULL);
     1089        int d3 = transform_distance_line_sq(c,d,a,NULL);
     1090        int d4 = transform_distance_line_sq(c,d,b,NULL);
     1091        if(!d1 || !d2 || !d3 || !d4) {
     1092                return 1;
     1093        }
     1094        return (geom_triangle_ccw(a,c,d) != geom_triangle_ccw(b,c,d)) && (geom_triangle_ccw(a,b,c) != geom_triangle_ccw(a,b,d));
     1095}
     1096
     1097int geom_poly_intersect(struct geom_poly_segment*a,struct geom_poly_segment*b)
     1098{
     1099        struct coord*curr_a = a->first;
     1100        struct coord*curr_b = b->first;
     1101        while (curr_a+1 <= a->last) {   //poly a segments
     1102                if(b->last!=b->first && geom_lineseg_intersect(curr_a,curr_a+1,b->last,b->first)) {
     1103                        return 1;
     1104                        break;
     1105                }
     1106                struct coord*curr_b = b->first;
     1107                while (curr_b+1 <= b->last) {   //poly b segments
     1108                        if(geom_lineseg_intersect(curr_a,curr_a+1,curr_b,curr_b+1)) {
     1109                                return 1;
     1110                                break;
     1111                        }
     1112                        curr_b++;
     1113                }
     1114                curr_a++;
     1115        }
     1116        if(b->last!=b->first && a->last!=a->first && geom_lineseg_intersect(a->last,a->first,b->last,b->first)) {
     1117                return 1;
     1118        }
     1119       
     1120        return 0;
     1121}
     1122
  • maptool/Makefile.am

     
    55endif
    66
    77AM_CPPFLAGS = @NAVIT_CFLAGS@ -I$(top_srcdir)/navit @ZLIB_CFLAGS@ @POSTGRESQL_CFLAGS@ -DMODULE=maptool
    8 libmaptool_la_SOURCES = boundaries.c buffer.c ch.c coastline.c geom.c itembin.c itembin_buffer.c misc.c osm.c osm_o5m.c osm_psql.c osm_protobuf.c osm_protobufdb.c osm_relations.c osm_xml.c sourcesink.c tempfile.c tile.c zip.c maptool.h generated-code/fileformat.pb-c.c generated-code/fileformat.pb-c.h generated-code/osmformat.pb-c.c generated-code/osmformat.pb-c.h google/protobuf-c/protobuf-c.c google/protobuf-c/protobuf-c.h google/protobuf-c/protobuf-c-private.h
     8libmaptool_la_SOURCES = boundaries.c buffer.c ch.c coastline.c geom.c itembin.c itembin_buffer.c misc.c osm.c osm_o5m.c osm_psql.c osm_protobuf.c osm_protobufdb.c osm_relations.c osm_xml.c sourcesink.c tempfile.c tile.c zip.c maptool.h generated-code/fileformat.pb-c.c generated-code/fileformat.pb-c.h generated-code/osmformat.pb-c.c generated-code/osmformat.pb-c.h google/protobuf-c/protobuf-c.c google/protobuf-c/protobuf-c.h google/protobuf-c/protobuf-c-private.h seidel/interface.h seidel/triangulate.h seidel/misc_seidel.c seidel/construct.c seidel/monotone.c  seidel/tri.c
    99maptool_SOURCES = maptool.c
    1010maptool_LDADD = libmaptool.la ../libnavit.la @NAVIT_LIBS@ @WORDEXP_LIBS@ @ZLIB_LIBS@ @POSTGRESQL_LIBS@ @CRYPTO_LIBS@ @INTLLIBS@ @LIBC_LIBS@
  • maptool/boundaries.c

     
    1919#include <stdio.h>
    2020#include <string.h>
    2121#include "maptool.h"
     22#include "types.h"
     23#include "seidel/triangulate.h"
    2224
    2325char *
    2426osm_tag_value(struct item_bin *ib, char *key)
     
    7072        while ((ib=read_item(boundaries))) {
    7173                char *member=NULL;
    7274                struct boundary *boundary=g_new0(struct boundary, 1);
     75                boundary->area_item_types = NULL;
     76                enum item_type* p_item_type;
     77                if((p_item_type = item_bin_get_attr(ib,attr_item_types,NULL))) {
     78                        int i=0;
     79                        do {
     80                                boundary->area_item_types = g_realloc(boundary->area_item_types, (i+1)*sizeof(enum item_type));
     81                                boundary->area_item_types[i] = *p_item_type;
     82                                i++;
     83                                boundary->area_item_types[i] = type_none;
     84                        } while(*p_item_type++!=type_none);
     85                }
    7386                char *admin_level=osm_tag_value(ib, "admin_level");
    7487                char *iso=osm_tag_value(ib, "ISO3166-1");
    7588                /* disable spain for now since it creates a too large index */
     
    192205        return g_list_prepend(list, boundary);
    193206}
    194207
     208//creates  inner attr for the first outer
     209struct attr*
     210create_triangle_attr(GList*grouped_segments)
     211{
     212        struct attr *ret = NULL;
     213        struct attr *triangle_attr = NULL;
     214        GList *l = grouped_segments;
     215        while (l) {
     216                struct geom_poly_segment *gs=l->data;
     217                if (coord_is_equal(*gs->first,*gs->last) && (gs->last-gs->first)>0) {
     218                        gs->last--;     
     219                }
     220                if(gs->type==geom_poly_segment_type_way_outer && (gs->last-gs->first+1)>2 ) {   //use the first outer
     221                        s32* polylen_buf = NULL;
     222                        s32  polylen = 0;
     223                        s32* vertices_buf = NULL;
     224                        GList *k = grouped_segments;
     225                        s32 prev_size = 0;
     226                        s32 coordsize_in_int=(sizeof(struct coord))/(sizeof(s32));
    195227
     228                        if(geom_polygon_is_clockwise(gs)) { // outer polygons should be CCW
     229                                geom_polygon_reverse(gs);
     230                        }
     231#if 0
     232printf("CRD:START\n");
     233fflush(stdout);
     234#endif
     235                        while(k) { //loop through inner polygons
     236                                struct geom_poly_segment *gs2=k->data;
     237                                //if actual inner poly inside actual outer, add inner to the list of to be added inner polygons
     238                                if (coord_is_equal(*gs2->first,*gs2->last) && (gs2->last-gs2->first)>0) {
     239                                        gs2->last--;   
     240                                }
     241
     242                                if(gs2->type==geom_poly_segment_type_way_inner && (gs2->last-gs2->first+1)>2 && geom_polygon_in_polygon(gs2, gs) && !geom_poly_intersect(gs2,gs)) {
     243                                        //make sure that outer is CCW, inner is CW
     244                                        if(!geom_polygon_is_clockwise(gs2)) { // inner polygons should be CW
     245                                                geom_polygon_reverse(gs2);
     246                                        }
     247                                        if(!polylen_buf) { //collect info about the outer
     248                                                polylen_buf = g_realloc(polylen_buf,sizeof(s32)*coordsize_in_int);
     249                                                polylen_buf[0] = gs->last-gs->first+1;
     250                                                ++polylen;
     251
     252#if 0
     253struct coord*c=gs->first;
     254do {
     255        printf("CRD:%d %d\n",c->x,c->y);
     256        c++;
     257} while(c<=gs->last);
     258printf("CRD:\n");
     259fflush(stdout);
     260#endif
     261                                                vertices_buf = g_realloc(vertices_buf,sizeof(s32)*(coordsize_in_int*(gs->last-gs->first+1)));
     262                                                prev_size += coordsize_in_int*(gs->last-gs->first+1);
     263                                                memcpy(vertices_buf,gs->first,sizeof(s32)*coordsize_in_int*(gs->last-gs->first+1));
     264                                        }
     265#if 0
     266struct coord*c=gs2->first;
     267do {
     268        printf("CRD:%d %d\n",c->x,c->y);
     269        c++;
     270} while(c<=gs2->last);
     271printf("CRD:\n");
     272fflush(stdout);
     273#endif
     274
     275                                        /*
     276                                        */
     277                                        polylen_buf = g_realloc(polylen_buf,sizeof(s32)*coordsize_in_int*(++polylen));
     278                                        polylen_buf[polylen-1] = gs2->last-gs2->first+1;
     279                                        vertices_buf = g_realloc(vertices_buf,sizeof(s32)*(prev_size+coordsize_in_int*(gs2->last-gs2->first+1)));
     280                                        memcpy(&vertices_buf[prev_size],gs2->first,sizeof(s32)*coordsize_in_int*(gs2->last-gs2->first+1));
     281                                        prev_size+=coordsize_in_int*(gs2->last-gs2->first+1);
     282                                }
     283                                k=g_list_next(k);
     284                        }
     285#if 0
     286printf("CRD:END\n");
     287fflush(stdout);
     288#endif
     289
     290struct coord*c = (struct coord*)vertices_buf;
     291struct coord*startcoord = (struct coord*)vertices_buf;
     292int duplicate=0;
     293while(c<(startcoord+prev_size/coordsize_in_int) && duplicate==0) {
     294        struct coord*c2 = c+1;
     295        while(c2<(startcoord+prev_size/coordsize_in_int)) {
     296                if(coord_is_equal(*c,*c2)) {
     297                        duplicate=1;
     298                        break;
     299                }
     300                c2++;
     301        }
     302        c++;
     303}
     304if( duplicate==0 && prev_size<SEGSIZE) {
     305                        s32 sum_size = prev_size;
     306                        s32 tri_num = (sum_size - 2) + 2*(polylen-1);
     307                        if(tri_num>0 && polylen>0) {
     308                                int ii=0;
     309                                s32 (*triangles)[3] = g_malloc0(sizeof(s32)*tri_num*3);
     310#if 0
     311printf("polylen %d\n",polylen);
     312for(ii=0;ii<polylen;++ii) {
     313        printf("polylen[%d] == %d\n",ii,polylen_buf[ii]);
     314}
     315fflush(stdout);
     316#endif
     317
     318                                int tri_ret = triangulate_polygon(polylen, polylen_buf, startcoord-1, triangles);
     319                                if(tri_ret==0) {
     320                                  //TODO free polylen_buf and vertices_buf when necessary
     321                                  g_free(polylen_buf);
     322                                  g_free(vertices_buf);
     323                                  g_free(triangles);
     324                                  return NULL;
     325                                }
     326
     327                                s32* pit=g_malloc(sizeof(s32)*(1 + 1 + 3*(tri_ret)) + sizeof(struct coord)*sum_size);
     328                                pit[0]=sum_size;
     329                                pit[1]=tri_ret;
     330                                memcpy(&pit[2],vertices_buf,sizeof(s32)*sum_size);
     331                                memcpy(&pit[2+sum_size],triangles,sizeof(s32)*3*tri_ret);
     332
     333                                triangle_attr = g_new0(struct attr,1);
     334                                s32* pi;
     335                                triangle_attr->type = attr_triangles;
     336                                triangle_attr->u.data = pit;
     337       
     338                                g_free(triangles);
     339                        }
     340
     341                        ret = triangle_attr;
     342
     343                        g_free(polylen_buf);
     344                        g_free(vertices_buf);
     345                        return ret;
     346                }
     347                g_free(polylen_buf);
     348                g_free(vertices_buf);
     349        }
     350
     351//TODO merge subsequent repeated points
     352
     353                l=g_list_next(l);
     354        }
     355        return NULL;
     356}
     357
     358/*
     359 * removes inner polygons intersecting the outer and other inners
     360 * expects a grouped polygon lists
     361 */
     362GList* remove_intersecting_segments(GList*segments)
     363{
     364        //find outer
     365        GList*curr=segments;
     366        GList*curr2=segments;
     367        struct geom_poly_segment*outer=NULL;
     368        int changed;
     369        int changed_g;
     370        while(curr) { //find outer
     371                struct geom_poly_segment *gs = curr->data;
     372                changed=0;
     373                if ( gs->last-gs->first<=0 || (gs->type==geom_poly_segment_type_way_outer && self_intersect_test_segment(gs)) ) {
     374                        curr=g_list_next(curr);
     375                        segments=g_list_remove(segments,gs);
     376                        changed=1;
     377                } else if(gs->type==geom_poly_segment_type_way_outer) {
     378                        outer=gs;
     379                }
     380                if(!changed) {
     381                        curr = g_list_next(curr);
     382                }
     383        }
     384
     385        curr=segments;
     386        while(outer && curr) {
     387                changed=0;
     388                struct geom_poly_segment *gs = curr->data;
     389                if(gs->type==geom_poly_segment_type_way_inner) {
     390                        if(geom_poly_intersect(outer,gs) || self_intersect_test_segment(gs) || geom_polygon_in_polygon(outer,gs)  ) {
     391                                curr=g_list_next(curr);
     392                                segments=g_list_remove(segments,gs);
     393                                changed=1;
     394                        }
     395                } else if(gs->type==geom_poly_segment_type_way_outer && gs!=outer) {
     396                }
     397                if(!changed) {
     398                        curr = g_list_next(curr);
     399                }
     400        }
     401       
     402        do {
     403                changed_g=0;
     404                curr=segments;
     405
     406        while(curr) {
     407                struct geom_poly_segment *gs = curr->data;
     408                curr2 = g_list_next(curr);
     409                while(curr2) {
     410                        struct geom_poly_segment *gs2 = curr2->data;
     411                        changed=0;
     412                        if(self_intersect_test_segment(gs2)) {
     413                                curr2=g_list_next(curr2);
     414                                segments=g_list_remove(segments,gs2);
     415                                changed=1;
     416                                changed_g=1;
     417                                continue;
     418                        }
     419                        if(gs->type==geom_poly_segment_type_way_inner && gs2->type==geom_poly_segment_type_way_inner) {
     420                                if( geom_polygon_in_polygon(gs2,gs) || /*currently we do not support nested inners*/
     421                                        geom_poly_intersect(gs,gs2)) {
     422
     423                                        curr2=g_list_next(curr2);
     424                                        segments=g_list_remove(segments,gs2);
     425                                        changed=1;
     426                                        changed_g=1;
     427                                }
     428                        }
     429                        if(!changed) {
     430                                curr2 = g_list_next(curr2);
     431                        }
     432                }
     433                curr = g_list_next(curr);
     434        }
     435
     436        } while(changed_g);
     437        return segments;
     438}
     439
    196440static GList *
    197441process_boundaries_finish(GList *boundaries_list)
    198442{
    199         GList *l,*sl,*l2,*ln;
     443        GList *l,*sl;
    200444        GList *ret=NULL;
    201445        l=boundaries_list;
    202446        while (l) {
    203447                struct boundary *boundary=l->data;
    204448                int first=1;
    205                 FILE *f=NULL,*fu=NULL;
     449                FILE *f=NULL,*fu=NULL,*ways_split=NULL;
     450                int to_write = 0;
     451                int have_inner = 0;
    206452                if (boundary->country) {
    207453                        char *name=g_strdup_printf("country_%s_poly",boundary->iso2);
    208454                        f=tempfile("",name,1);
    209455                        g_free(name);
    210456                }
     457                if (boundary->area_item_types) {
     458                        ways_split=tempfile("","ways_split",2);
     459                        struct item_bin *ib=item_bin_relation_area;
     460                        item_bin_init(ib, boundary->area_item_types[0]);
     461                }
    211462                boundary->sorted_segments=geom_poly_segments_sort(boundary->segments, geom_poly_segment_type_way_right_side);
    212463                sl=boundary->sorted_segments;
    213464                while (sl) {
    214465                        struct geom_poly_segment *gs=sl->data;
     466                        geom_polygon_unique(gs);
     467
    215468                        struct coord *c=gs->first;
    216469                        while (c <= gs->last) {
    217470                                if (first) {
     
    228481                                item_bin_add_coord(ib, gs->first, gs->last-gs->first+1);
    229482                                item_bin_write(ib, f);
    230483                        }
     484                        if (ways_split && gs->type==geom_poly_segment_type_way_outer) {
     485                                struct item_bin *ib=item_bin_relation_area;
     486                                item_bin_add_coord(ib, gs->first, gs->last-gs->first+1);
     487                                to_write = 1;
     488                        }
     489                        if (ways_split && gs->type==geom_poly_segment_type_way_inner) {
     490                                have_inner = 1;
     491                        }
    231492                        if (boundary->country) {
    232493                                if (!coord_is_equal(*gs->first,*gs->last)) {
    233494                                        if (!fu) {
     
    246507                        }
    247508                        sl=g_list_next(sl);
    248509                }       
     510                if (ways_split && to_write ) {
     511                        struct item_bin *ib=item_bin_relation_area;
     512                        if(0 == self_intersect_test(ib)) { //one outer
     513                                int i=0;
     514                                        if(have_inner) {
     515                                                GList*grouped_segments =
     516                                                        geom_poly_segments_group(boundary->sorted_segments, grouped_segments);
     517                                                grouped_segments = remove_intersecting_segments(grouped_segments);
     518#if 1
     519printf("OSMID:%d\n",item_bin_get_relationid(boundary->ib));
     520fflush(stdout);
     521#endif
     522                                                struct attr*triangle_attr = create_triangle_attr(grouped_segments);
     523
     524                                                if(triangle_attr && triangle_attr->u.data) {
     525                                                        item_bin_add_attr(ib,triangle_attr);
     526                                                }
     527                                                if(triangle_attr) {
     528                                                        g_free(triangle_attr->u.data);
     529                                                        g_free(triangle_attr);
     530                                                }
     531                                        }
     532
     533                                i = 0;
     534                                while(boundary->area_item_types[i]!=type_none) { //write all matches
     535                                        ib->type = boundary->area_item_types[i++];
     536                                        item_bin_write(ib, ways_split);
     537                                };
     538                        } else { /*maybe disjunct outer polygons*/
     539                                GList*grouped_segments =
     540                                        geom_poly_segments_group(boundary->sorted_segments, grouped_segments);
     541                                grouped_segments = remove_intersecting_segments(grouped_segments);
     542
     543                                GList*l = grouped_segments;
     544                                int self_intersecting_grouped = 0;
     545                                int i=0;
     546                                int changed=0;
     547                                while (l) {
     548                                        int changed=0;
     549                                        struct geom_poly_segment *gs=l->data;
     550                                                if(gs->type==geom_poly_segment_type_way_outer) {
     551                                                struct item_bin *ib=item_bin_relation_area;
     552                                                item_bin_init(ib,boundary->area_item_types[0]);
     553
     554                                                if(gs->last-gs->first+1 > 0) {
     555                                                                item_bin_add_coord(ib, gs->first, gs->last-gs->first+1);
     556                                                }
     557                                                if( 0!=self_intersect_test(ib)) {
     558                                                        self_intersecting_grouped = 1;
     559                                                } else {
     560#if 0
     561printf("OSMID2:%d   boundary->area_item_types[0]=%X\n",item_bin_get_relationid(boundary->ib),boundary->area_item_types[0]);
     562fflush(stdout);
     563#endif
     564                                                        struct attr *triangle_attr = create_triangle_attr(grouped_segments);
     565
     566                                                        if(triangle_attr  && triangle_attr->u.data) {
     567                                                                int i=0;
     568                                                                item_bin_add_attr(ib,triangle_attr);
     569                                                                g_free(triangle_attr->u.data);
     570                                                        }
     571                                                        i=0;
     572                                                        while(boundary->area_item_types[i]!=type_none) { //write all match
     573                                                                ib->type = boundary->area_item_types[i];
     574                                                                item_bin_write(ib, ways_split);
     575                                                                ++i;
     576                                                        };
     577                                                }
     578                                                //remove first outer
     579                                                changed=1;
     580                                                l=g_list_next(l);
     581                                                grouped_segments = g_list_remove(grouped_segments,gs);
     582                                        }
     583                                        if(!changed) {
     584                                                l=g_list_next(l);
     585                                        }
     586                                }
     587
     588                                if(self_intersecting_grouped) {
     589                                        osm_warning("relation",item_bin_get_relationid(boundary->ib),0,"Self intersecting relation area\n");
     590                                }
     591                                g_list_free(grouped_segments);
     592                        }
     593                }
     594
    249595                ret=process_boundaries_insert(ret, boundary);
    250596                l=g_list_next(l);
    251597                if (f)
    252598                        fclose(f);
     599                if (ways_split)
     600                        fclose(ways_split);
    253601                if (fu) {
    254602                        if (boundary->country)
    255603                                osm_warning("relation",item_bin_get_relationid(boundary->ib),0,"Broken country polygon '%s'\n",boundary->iso2);
    256604                        fclose(fu);
    257605                }
    258                
     606        g_free(boundary->area_item_types);
    259607        }
    260608#if 0
    261609        printf("hierarchy\n");
  • maptool/CMakeLists.txt

     
    22if(BUILD_MAPTOOL)
    33   add_definitions( -DMODULE=maptool ${NAVIT_COMPILE_FLAGS})
    44   include_directories(${CMAKE_CURRENT_SOURCE_DIR})
    5    SET(MAPTOOL_SOURCE boundaries.c buffer.c ch.c coastline.c geom.c itembin.c itembin_buffer.c misc.c osm.c osm_o5m.c osm_relations.c sourcesink.c tempfile.c tile.c zip.c osm_xml.c)
     5   SET(MAPTOOL_SOURCE boundaries.c buffer.c ch.c coastline.c geom.c itembin.c itembin_buffer.c misc.c osm.c osm_o5m.c osm_relations.c sourcesink.c tempfile.c tile.c zip.c osm_xml.c seidel/misc_seidel.c seidel/construct.c seidel/monotone.c  seidel/tri.c)
    66   if(NOT MSVC)
    77        SET(MAPTOOL_SOURCE ${MAPTOOL_SOURCE} osm_protobuf.c osm_protobufdb.c generated-code/fileformat.pb-c.c generated-code/osmformat.pb-c.c google/protobuf-c/protobuf-c.c)
    88   endif(NOT MSVC)
  • maptool/osm.c

     
    672672        "w      natural=scree           poly_scree\n"
    673673        "w      natural=scrub           poly_scrub\n"
    674674        "w      natural=water           poly_water\n"
     675        "w      natural=wetland         poly_wetland\n"
    675676        "w      natural=wood            poly_wood\n"
    676677        "w      piste:type=downhill,piste:difficulty=advanced           piste_downhill_advanced\n"
    677678        "w      piste:type=downhill,piste:difficulty=easy               piste_downhill_easy\n"
     
    709710        "w      waterway=canal          water_canal\n"
    710711        "w      waterway=drain          water_drain\n"
    711712        "w      waterway=river          water_river\n"
     713        "w      water=river             water_river\n"
    712714        "w      waterway=riverbank      poly_water\n"
    713715        "w      waterway=stream         water_stream\n"
    714716        "w      barrier=ditch   ditch\n"
     
    942944        char buffer[BUFFER_SIZE*2+2];
    943945        if (in_relation) {
    944946                relation_add_tag(k,v);
     947
     948                //for relation areas we don't set flags
     949                strcpy(buffer,"*=*");
     950                if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer)))
     951                        attr_present[idx]=1;
     952
     953                sprintf(buffer,"%s=*", k);
     954                if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer)))
     955                        attr_present[idx]=2;
     956
     957                sprintf(buffer,"*=%s", v);
     958                if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer)))
     959                        attr_present[idx]=2;
     960
     961                sprintf(buffer,"%s=%s", k, v);
     962                if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer)))
     963                        attr_present[idx]=4;
     964
    945965                return;
    946966        }
    947967        if (! strcmp(k,"ele"))
     
    14821502}
    14831503
    14841504
     1505static int
     1506attr_longest_match(struct attr_mapping **mapping, int mapping_count, enum item_type *types, int types_count)
     1507{
     1508        int i,j,longest=0,ret=0,sum,val;
     1509        struct attr_mapping *curr;
     1510        for (i = 0 ; i < mapping_count ; i++) {
     1511                sum=0;
     1512                curr=mapping[i];
     1513                for (j = 0 ; j < curr->attr_present_idx_count ; j++) {
     1514                        val=attr_present[curr->attr_present_idx[j]];
     1515                        if (val)
     1516                                sum+=val;
     1517                        else {
     1518                                sum=-1;
     1519                                break;
     1520                        }
     1521                }
     1522                if (sum > longest) {
     1523                        longest=sum;
     1524                        ret=0;
     1525                }
     1526                if (sum > 0 && sum == longest && ret < types_count)
     1527                        types[ret++]=curr->type;
     1528        }
     1529        return ret;
     1530}
     1531
     1532static void
     1533attr_longest_match_clear(void)
     1534{
     1535        memset(attr_present, 0, sizeof(*attr_present)*attr_present_count);
     1536}
     1537
    14851538void
    14861539osm_end_relation(struct maptool_osm *osm)
    14871540{
     1541        int count,itcount=0;
     1542        enum item_type types[10];
     1543        types[0] = type_none;
     1544
     1545        count=attr_longest_match(attr_mapping_way, attr_mapping_way_count, types, sizeof(types)/sizeof(enum item_type));
     1546        if(count>0) {
     1547                int ii;
     1548
     1549                struct attr attr;
     1550                attr.type = attr_item_types;
     1551                attr.u.item_types = NULL;
     1552
     1553                for(ii=0;ii<count;ii++) {
     1554                        if(item_type_is_area(types[ii])) {
     1555                                boundary = 1;
     1556                                attr.u.item_types = g_realloc(attr.u.item_types, (itcount+2)*sizeof(enum item_type));
     1557                                attr.u.item_types[itcount] = types[ii];
     1558                                attr.u.item_types[itcount+1] = type_none;
     1559                                ++itcount;
     1560                        }
     1561                }
     1562                if(itcount>0) {
     1563                        item_bin_add_attr(item_bin, &attr);
     1564                }
     1565        }
     1566        attr_longest_match_clear();
     1567
    14881568        in_relation=0;
    14891569        if ((!strcmp(relation_type, "multipolygon") || !strcmp(relation_type, "boundary")) && boundary) {
    14901570#if 0
     
    15531633        }
    15541634}
    15551635
    1556 
    1557 static int
    1558 attr_longest_match(struct attr_mapping **mapping, int mapping_count, enum item_type *types, int types_count)
    1559 {
    1560         int i,j,longest=0,ret=0,sum,val;
    1561         struct attr_mapping *curr;
    1562         for (i = 0 ; i < mapping_count ; i++) {
    1563                 sum=0;
    1564                 curr=mapping[i];
    1565                 for (j = 0 ; j < curr->attr_present_idx_count ; j++) {
    1566                         val=attr_present[curr->attr_present_idx[j]];
    1567                         if (val)
    1568                                 sum+=val;
    1569                         else {
    1570                                 sum=-1;
    1571                                 break;
    1572                         }
    1573                 }
    1574                 if (sum > longest) {
    1575                         longest=sum;
    1576                         ret=0;
    1577                 }
    1578                 if (sum > 0 && sum == longest && ret < types_count)
    1579                         types[ret++]=curr->type;
    1580         }
    1581         return ret;
    1582 }
    1583 
    1584 static void
    1585 attr_longest_match_clear(void)
    1586 {
    1587         memset(attr_present, 0, sizeof(*attr_present)*attr_present_count);
    1588 }
    1589 
    15901636void
    15911637osm_end_way(struct maptool_osm *osm)
    15921638{
     
    25622608        processed_nodes=processed_nodes_out=processed_ways=processed_relations=processed_tiles=0;
    25632609        sig_alrm(0);
    25642610        while ((ib=read_item(in))) {
     2611                int self_intersect;
    25652612#if 0
    25662613                fprintf(stderr,"type 0x%x len %d clen %d\n", ib->type, ib->len, ib->clen);
    25672614#endif
     
    25772624                                if (ni) {
    25782625                                        c[i]=ni->c;
    25792626                                        if (ni->ref_way > 1 && i != 0 && i != ccount-1 && i != last && item_get_default_flags(ib->type)) {
    2580                                                 write_item_part(out, out_index, out_graph, ib, last, i, &last_id);
     2627                                                if(!(0 != self_intersect_test(ib) && item_type_is_area(ib->type)))
     2628                                                        write_item_part(out, out_index, out_graph, ib, last, i, &last_id);
    25812629                                                last=i;
    25822630                                        }
    25832631                                } else if (final) {
     
    25922640                                }
    25932641                        }
    25942642                }
     2643                self_intersect = self_intersect_test(ib);
     2644                if(self_intersect!=0 && (/*ib->type == type_water_line ||*/ item_type_is_area(ib->type)) ) {
     2645                        osm_warning("way", item_bin_get_wayid(ib), 0, "Self intersecting area\n");
     2646                }
    25952647                if (ccount) {
    2596                         write_item_part(out, out_index, out_graph, ib, last, ccount-1, &last_id);
     2648                        if(!(0 != self_intersect && item_type_is_area(ib->type)))
     2649                                write_item_part(out, out_index, out_graph, ib, last, ccount-1, &last_id);
    25972650                        if (final && ib->type == type_water_line && out_coastline) {
    2598                                 write_item_part(out_coastline, NULL, NULL, ib, last, ccount-1, NULL);
     2651                                //if(0 == self_intersect)
     2652                                        write_item_part(out_coastline, NULL, NULL, ib, last, ccount-1, NULL);
    25992653                        }
    26002654                }
    26012655        }
  • maptool/maptool.c

     
    745745        }
    746746}
    747747
    748 
    749748int main(int argc, char **argv)
    750749{
    751750#if 0
     
    821820                return 0;
    822821        }
    823822        phase=0;
    824 
     823/**/
    825824    // input from an OSM file
    826825        if (p.input == 0) {
    827826                if (start_phase(&p, "collecting data")) {
     
    857856        if (start_phase(&p,"generating coastlines")) {
    858857                osm_process_coastlines(&p, suffix);
    859858        }
     859/**/
    860860        if (start_phase(&p,"assinging towns to countries")) {
    861861                FILE *towns=tempfile(suffix,"towns",0),*boundaries=NULL,*ways=NULL;
    862862                if (towns) {
  • attr_def.h

     
    430430ATTR(announcement)
    431431ATTR(cursor)
    432432ATTR(config)
     433ATTR(triangles)
    433434ATTR2(0x0008ffff,type_object_end)
    434435ATTR2(0x00090000,type_coord_begin)
    435436ATTR2(0x0009ffff,type_coord_end)
  • attr.c

     
    492492int
    493493attr_data_size(struct attr *attr)
    494494{
     495        if (attr->type == attr_triangles) {
     496                int*pi=(u32*)(attr->u.data);
     497                return sizeof(u32)*(2 + 3*pi[1] + pi[0]);
     498        }
    495499        if (attr->type >= attr_type_string_begin && attr->type <= attr_type_string_end)
    496500                return strlen(attr->u.str)+1;
    497501        if (attr->type >= attr_type_int_begin && attr->type <= attr_type_int_end)
  • graphics.c

     
    5050#include "callback.h"
    5151#include "file.h"
    5252#include "event.h"
     53#include "types.h"
    5354
    5455
    5556//##############################################################################################################
     
    452453*/
    453454void graphics_gc_destroy(struct graphics_gc *gc)
    454455{
    455         if (!gc)
    456             return;
    457456        gc->meth.gc_destroy(gc->priv);
    458457        g_free(gc);
    459458}
     
    860859        char *label;
    861860        int z_order;
    862861        int count;
     862        struct attr *triangles;/*stores triangle data for polygons that need triangulation (polygons with holes)*/
    863863        struct coord c[0];
    864864};
    865865
     
    876876                struct displayitem *di=dl->hash_entries[i].di;
    877877                while (di) {
    878878                        struct displayitem *next=di->next;
     879                        g_free(di->triangles);
    879880                        g_free(di);
    880881                        di=next;
    881882                }
     
    889890 * @returns <>
    890891 * @author Martin Schaller (04/2008)
    891892*/
    892 static void display_add(struct hash_entry *entry, struct item *item, int count, struct coord *c, char **label, int label_count)
     893static void display_add(struct hash_entry *entry, struct item *item, int count, struct coord *c, char **label, int label_count, struct attr*poly_attr)
    893894{
    894895        struct displayitem *di;
    895896        int len,i;
     
    910911        p+=sizeof(*di)+count*sizeof(*c);
    911912        di->item=*item;
    912913        di->z_order=0;
     914        di->triangles = poly_attr;
    913915        if (label && label_count) {
    914916                di->label=p;
    915917                for (i = 0 ; i < label_count ; i++) {
     
    17271729        char *path;
    17281730
    17291731        while (di) {
     1732
    17301733        int i,count=di->count,mindist=dc->mindist;
    17311734
    17321735        di->z_order=++(gra->current_z_order);
     
    17411744        if (dc->type == type_poly_water_tiled)
    17421745                mindist=0;
    17431746        if (dc->e->type == element_polyline)
    1744                 count=transform(dc->trans, dc->pro, di->c, pa, count, mindist, e->u.polyline.width, width);
     1747                count=transform(dc->trans, dc->pro, di->c, pa, count, 0/*mindist*/, e->u.polyline.width, width);
    17451748        else
    1746                 count=transform(dc->trans, dc->pro, di->c, pa, count, mindist, 0, NULL);
     1749                count=transform(dc->trans, dc->pro, di->c, pa, count, 0/*mindist*/, 0, NULL);
    17471750        switch (e->type) {
    17481751        case element_polygon:
    1749                 graphics_draw_polygon_clipped(gra, gc, pa, count);
     1752                if(di->triangles) {
     1753                        u32*pi=di->triangles->u.data;
     1754                        int ii;
     1755#if 0
     1756int jj;
     1757for (jj=0;jj<pi[1]*3+pi[0];jj++) {
     1758        printf("pi[%d]=%d\n",jj,pi[jj]);
     1759}
     1760#endif
     1761        struct coord* vertices=&pi[2];
     1762        --vertices;
     1763        int coord_size_in_int = sizeof(struct coord)/sizeof(int);
     1764
     1765        for(ii=0;ii<pi[1];++ii) {
     1766                struct coord ca[3];
     1767                struct point pa[3];
     1768                ca[0].x = vertices[ pi[ 2 + pi[0] + ii*3 + 0 ] ].x;
     1769                ca[0].y = vertices[ pi[ 2 + pi[0] + ii*3 + 0 ] ].y;
     1770                ca[1].x = vertices[ pi[ 2 + pi[0] + ii*3 + 1 ] ].x;
     1771                ca[1].y = vertices[ pi[ 2 + pi[0] + ii*3 + 1 ] ].y;
     1772                ca[2].x = vertices[ pi[ 2 + pi[0] + ii*3 + 2 ] ].x;
     1773                ca[2].y = vertices[ pi[ 2 + pi[0] + ii*3 + 2 ] ].y;
     1774                count=transform(dc->trans, dc->pro, ca, pa, 3, 0, 0, NULL);
     1775                //TODO use some caching mechanism to avoid unnecessary recaculations
     1776                /*
     1777                if(
     1778                        (pa[0].x<gra->r.lu.x || gra->r.rl.x<pa[0].x || pa[0].y<gra->r.lu.y || gra->r.rl.y<pa[0].y) &&
     1779                        (pa[1].x<gra->r.lu.x || gra->r.rl.x<pa[1].x || pa[1].y<gra->r.lu.y || gra->r.rl.y<pa[1].y) &&
     1780                        (pa[2].x<gra->r.lu.x || gra->r.rl.x<pa[2].x || pa[2].y<gra->r.lu.y || gra->r.rl.y<pa[2].y)
     1781               
     1782                ) {
     1783                        //continue;
     1784                }*/
     1785                graphics_draw_polygon_clipped(gra, gc, &pa[0], 3);
     1786        }
     1787
     1788                } else {
     1789                        graphics_draw_polygon_clipped(gra, gc, pa, count);
     1790                }
    17501791                break;
    17511792        case element_polyline:
    17521793                {
     
    20892130                }
    20902131                if (displaylist->mr) {
    20912132                        while ((item=map_rect_get_item(displaylist->mr))) {
     2133
     2134
     2135
     2136
     2137
     2138
     2139
    20922140                                int label_count=0;
    20932141                                char *labels[2];
    20942142                                struct hash_entry *entry;
     
    21102158                                        dbg(0,"point count overflow %d for %s "ITEM_ID_FMT"\n", count,item_to_name(item->type),ITEM_ID_ARGS(*item));
    21112159                                        displaylist->dc.maxlen=max*2;
    21122160                                }
     2161
     2162                struct attr* poly_attr = g_new0(struct attr,1);
     2163
     2164                if(item_attr_get(item,attr_triangles,poly_attr)) {
     2165                        int *pi = poly_attr->u.data;
     2166                        int *pi2;
     2167                        int ii;
     2168                        pi2 = g_malloc(attr_data_size(poly_attr));
     2169                        memcpy(pi2,pi,attr_data_size(poly_attr));
     2170                        poly_attr->u.data = pi2;
     2171
     2172                } else {
     2173                        g_free(poly_attr);
     2174                        poly_attr = NULL;
     2175                }
     2176
    21132177                                if (item_is_custom_poi(*item)) {
    21142178                                        if (item_attr_get(item, attr_icon_src, &attr2))
    21152179                                                labels[1]=map_convert_string(displaylist->m, attr2.u.str);
     
    21282192                                        labels[0]=NULL;
    21292193                                if (displaylist->conv && label_count) {
    21302194                                        labels[0]=map_convert_string(displaylist->m, labels[0]);
    2131                                         display_add(entry, item, count, ca, labels, label_count);
     2195                                        display_add(entry, item, count, ca, labels, label_count ,poly_attr);
    21322196                                        map_convert_free(labels[0]);
    21332197                                } else
    2134                                         display_add(entry, item, count, ca, labels, label_count);
     2198                                        display_add(entry, item, count, ca, labels, label_count ,poly_attr);
    21352199                                if (labels[1])
    21362200                                        map_convert_free(labels[1]);
    21372201                                workload++;