Index: maptool/maptool.h =================================================================== --- maptool/maptool.h (revision 5005) +++ maptool/maptool.h (working copy) @@ -121,6 +121,7 @@ struct boundary { struct item_bin *ib; struct country_table *country; + enum item_type* area_item_types; char *iso2; GList *segments,*sorted_segments; GList *children; @@ -187,8 +188,14 @@ GList *geom_poly_segments_sort(GList *in, enum geom_poly_segment_type type); struct geom_poly_segment *item_bin_to_poly_segment(struct item_bin *ib, int type); int geom_poly_segments_point_inside(GList *in, struct coord *c); +GList* geom_poly_segments_group(GList *in, GList *out); void clip_line(struct item_bin *ib, struct rect *r, struct tile_parameter *param, struct item_bin_sink *out); void clip_polygon(struct item_bin *ib, struct rect *r, struct tile_parameter *param, struct item_bin_sink *out); +int self_intersect_test(struct item_bin*ib); +int geom_point_in_polygon(struct coord*c, struct geom_poly_segment*poly); +int geom_polygon_in_polygon(struct geom_poly_segment*poly1, struct geom_poly_segment*poly2); +int geom_polygon_is_clockwise(struct geom_poly_segment*poly); +void geom_polygon_reverse(struct geom_poly_segment*poly); /* itembin.c */ @@ -240,6 +247,7 @@ extern struct buffer node_buffer; extern int processed_nodes, processed_nodes_out, processed_ways, processed_relations, processed_tiles; extern struct item_bin *item_bin; +extern struct item_bin *item_bin_relation_area; extern int bytes_read; extern int overlap; extern int unknown_country; Index: maptool/itembin_buffer.c =================================================================== --- maptool/itembin_buffer.c (revision 5005) +++ maptool/itembin_buffer.c (working copy) @@ -24,6 +24,8 @@ static char buffer[800000]; struct item_bin *item_bin=(struct item_bin *)(void *)buffer; +static char buffer_relation_area[8000000]; +struct item_bin *item_bin_relation_area=(struct item_bin *)(void *)buffer_relation_area; static struct node_item *node_item=(struct node_item *)(void *)buffer; struct node_item * Index: maptool/seidel/tri.c =================================================================== --- maptool/seidel/tri.c (revision 0) +++ maptool/seidel/tri.c (revision 0) @@ -0,0 +1,166 @@ +#include "triangulate.h" +#include + + +static int initialise(n) + int n; +{ + register int i; + + for (i = 1; i <= n; i++) + seg[i].is_inserted = FALSE; + + generate_random_ordering(n); + + return 0; +} + +#ifdef STANDALONE + +int main(argc, argv) + int argc; + char *argv[]; +{ + int n, nmonpoly, genus; + int op[SEGSIZE][3], i, ntriangles; + + if ((argc < 2) || ((n = read_segments(argv[1], &genus)) < 0)) + { + fprintf(stderr, "usage: triangulate \n"); + exit(1); + } + + initialise(n); + construct_trapezoids(n); + nmonpoly = monotonate_trapezoids(n); + ntriangles = triangulate_monotone_polygons(n, nmonpoly, op); + + for (i = 0; i < ntriangles; i++) + printf("triangle #%d: %d %d %d\n", i, + op[i][0], op[i][1], op[i][2]); + + return 0; +} + + +#else /* Not standalone. Use this as an interface routine */ + + +/* Input specified as contours. + * Outer contour must be anti-clockwise. + * All inner contours must be clockwise. + * + * Every contour is specified by giving all its points in order. No + * point shoud be repeated. i.e. if the outer contour is a square, + * only the four distinct endpoints shopudl be specified in order. + * + * ncontours: #contours + * cntr: An array describing the number of points in each + * contour. Thus, cntr[i] = #points in the i'th contour. + * vertices: Input array of vertices. Vertices for each contour + * immediately follow those for previous one. Array location + * vertices[0] must NOT be used (i.e. i/p starts from + * vertices[1] instead. The output triangles are + * specified w.r.t. the indices of these vertices. + * triangles: Output array to hold triangles. + * + * Enough space must be allocated for all the arrays before calling + * this routine + */ + + +int triangulate_polygon(ncontours, cntr, vertices, triangles) + s32 ncontours; + s32 cntr[]; + s32 (*vertices)[2]; + s32 (*triangles)[3]; +{ + register int i; + int nmonpoly, ccount, npoints, genus; + int n; + + memset((void *)seg, 0, sizeof(seg)); + ccount = 0; + i = 1; + while (ccount < ncontours) + { + int j; + int first, last; + + npoints = cntr[ccount]; + first = i; + last = first + npoints - 1; + for (j = 0; j < npoints; j++, i++) + { + seg[i].v0.x = (double)(vertices[i][0]); + seg[i].v0.y = (double)vertices[i][1]; + + if (i == last) + { + seg[i].next = first; + seg[i].prev = i-1; + seg[i-1].v1 = seg[i].v0; + } + else if (i == first) + { + seg[i].next = i+1; + seg[i].prev = last; + seg[last].v1 = seg[i].v0; + } + else + { + seg[i].prev = i-1; + seg[i].next = i+1; + seg[i-1].v1 = seg[i].v0; + } + + seg[i].is_inserted = FALSE; + } + + ccount++; + } + + genus = ncontours - 1; + n = i-1; + + initialise(n); + if(construct_trapezoids(n)) { + return 0; + } + nmonpoly = monotonate_trapezoids(n); + return triangulate_monotone_polygons(n, nmonpoly, triangles); +} + + +/* This function returns TRUE or FALSE depending upon whether the + * vertex is inside the polygon or not. The polygon must already have + * been triangulated before this routine is called. + * This routine will always detect all the points belonging to the + * set (polygon-area - polygon-boundary). The return value for points + * on the boundary is not consistent!!! + */ + +int is_point_inside_polygon(vertex) + double vertex[2]; +{ + point_t v; + int trnum, rseg; + trap_t *t; + + v.x = vertex[0]; + v.y = vertex[1]; + + trnum = locate_endpoint(&v, &v, 1); + t = &tr[trnum]; + + if (t->state == ST_INVALID) + return FALSE; + + if ((t->lseg <= 0) || (t->rseg <= 0)) + return FALSE; + rseg = t->rseg; + return _greater_than_equal_to(&seg[rseg].v1, &seg[rseg].v0); +} + + +#endif /* STANDALONE */ Index: maptool/seidel/monotone.c =================================================================== --- maptool/seidel/monotone.c (revision 0) +++ maptool/seidel/monotone.c (revision 0) @@ -0,0 +1,738 @@ +#include "triangulate.h" +#include + +#define CROSS_SINE(v0, v1) ((v0).x * (v1).y - (v1).x * (v0).y) +#define LENGTH(v0) (sqrt((v0).x * (v0).x + (v0).y * (v0).y)) + + +static monchain_t mchain[TRSIZE]; /* Table to hold all the monotone */ + /* polygons . Each monotone polygon */ + /* is a circularly linked list */ + +static vertexchain_t vert[SEGSIZE]; /* chain init. information. This */ + /* is used to decide which */ + /* monotone polygon to split if */ + /* there are several other */ + /* polygons touching at the same */ + /* vertex */ + +static int mon[SEGSIZE]; /* contains position of any vertex in */ + /* the monotone chain for the polygon */ +static int visited[TRSIZE]; +static int chain_idx, op_idx, mon_idx; + + +static int triangulate_single_polygon(int, int, int, int (*)[3]); +static int traverse_polygon(int, int, int, int); + +/* Function returns TRUE if the trapezoid lies inside the polygon */ +static int inside_polygon(t) + trap_t *t; +{ + int rseg = t->rseg; + + if (t->state == ST_INVALID) + return 0; + + if ((t->lseg <= 0) || (t->rseg <= 0)) + return 0; + + if (((t->u0 <= 0) && (t->u1 <= 0)) || + ((t->d0 <= 0) && (t->d1 <= 0))) /* triangle */ + return (_greater_than(&seg[rseg].v1, &seg[rseg].v0)); + + return 0; +} + + +/* return a new mon structure from the table */ +static int newmon() +{ + return ++mon_idx; +} + + +/* return a new chain element from the table */ +static int new_chain_element() +{ + return ++chain_idx; +} + + +static double get_angle(vp0, vpnext, vp1) + point_t *vp0; + point_t *vpnext; + point_t *vp1; +{ + point_t v0, v1; + + v0.x = vpnext->x - vp0->x; + v0.y = vpnext->y - vp0->y; + + v1.x = vp1->x - vp0->x; + v1.y = vp1->y - vp0->y; + + if (CROSS_SINE(v0, v1) >= 0) /* sine is positive */ + return DOT(v0, v1)/LENGTH(v0)/LENGTH(v1); + else + return (-1.0 * DOT(v0, v1)/LENGTH(v0)/LENGTH(v1) - 2); +} + + +/* (v0, v1) is the new diagonal to be added to the polygon. Find which */ +/* chain to use and return the positions of v0 and v1 in p and q */ +static int get_vertex_positions(v0, v1, ip, iq) + int v0; + int v1; + int *ip; + int *iq; +{ + vertexchain_t *vp0, *vp1; + register int i; + double angle, temp; + int tp, tq; + + vp0 = &vert[v0]; + vp1 = &vert[v1]; + + /* p is identified as follows. Scan from (v0, v1) rightwards till */ + /* you hit the first segment starting from v0. That chain is the */ + /* chain of our interest */ + + angle = -4.0; + for (i = 0; i < 4; i++) + { + if (vp0->vnext[i] <= 0) + continue; + if ((temp = get_angle(&vp0->pt, &(vert[vp0->vnext[i]].pt), + &vp1->pt)) > angle) + { + angle = temp; + tp = i; + } + } + + *ip = tp; + + /* Do similar actions for q */ + + angle = -4.0; + for (i = 0; i < 4; i++) + { + if (vp1->vnext[i] <= 0) + continue; + if ((temp = get_angle(&vp1->pt, &(vert[vp1->vnext[i]].pt), + &vp0->pt)) > angle) + { + angle = temp; + tq = i; + } + } + + *iq = tq; + + return 0; +} + + +/* v0 and v1 are specified in anti-clockwise order with respect to + * the current monotone polygon mcur. Split the current polygon into + * two polygons using the diagonal (v0, v1) + */ +static int make_new_monotone_poly(mcur, v0, v1) + int mcur; + int v0; + int v1; +{ + int p, q, ip, iq; + int mnew = newmon(); + int i, j, nf0, nf1; + vertexchain_t *vp0, *vp1; + + vp0 = &vert[v0]; + vp1 = &vert[v1]; + + get_vertex_positions(v0, v1, &ip, &iq); + + p = vp0->vpos[ip]; + q = vp1->vpos[iq]; + + /* At this stage, we have got the positions of v0 and v1 in the */ + /* desired chain. Now modify the linked lists */ + + i = new_chain_element(); /* for the new list */ + j = new_chain_element(); + + mchain[i].vnum = v0; + mchain[j].vnum = v1; + + mchain[i].next = mchain[p].next; + mchain[mchain[p].next].prev = i; + mchain[i].prev = j; + mchain[j].next = i; + mchain[j].prev = mchain[q].prev; + mchain[mchain[q].prev].next = j; + + mchain[p].next = q; + mchain[q].prev = p; + + nf0 = vp0->nextfree; + nf1 = vp1->nextfree; + + vp0->vnext[ip] = v1; + + vp0->vpos[nf0] = i; + vp0->vnext[nf0] = mchain[mchain[i].next].vnum; + vp1->vpos[nf1] = j; + vp1->vnext[nf1] = v0; + + vp0->nextfree++; + vp1->nextfree++; + +#ifdef DEBUG + fprintf(stderr, "make_poly: mcur = %d, (v0, v1) = (%d, %d)\n", + mcur, v0, v1); + fprintf(stderr, "next posns = (p, q) = (%d, %d)\n", p, q); +#endif + + mon[mcur] = p; + mon[mnew] = i; + return mnew; +} + +/* Main routine to get monotone polygons from the trapezoidation of + * the polygon. + */ + +int monotonate_trapezoids(n) + int n; +{ + register int i; + int tr_start; + + memset((void *)vert, 0, sizeof(vert)); + memset((void *)visited, 0, sizeof(visited)); + memset((void *)mchain, 0, sizeof(mchain)); + memset((void *)mon, 0, sizeof(mon)); + + /* First locate a trapezoid which lies inside the polygon */ + /* and which is triangular */ + for (i = 0; i < TRSIZE; i++) + if (inside_polygon(&tr[i])) + break; + tr_start = i; + + /* Initialise the mon data-structure and start spanning all the */ + /* trapezoids within the polygon */ + +#if 0 + for (i = 1; i <= n; i++) + { + mchain[i].prev = i - 1; + mchain[i].next = i + 1; + mchain[i].vnum = i; + vert[i].pt = seg[i].v0; + vert[i].vnext[0] = i + 1; /* next vertex */ + vert[i].vpos[0] = i; /* locn. of next vertex */ + vert[i].nextfree = 1; + } + mchain[1].prev = n; + mchain[n].next = 1; + vert[n].vnext[0] = 1; + vert[n].vpos[0] = n; + chain_idx = n; + mon_idx = 0; + mon[0] = 1; /* position of any vertex in the first */ + /* chain */ + +#else + + for (i = 1; i <= n; i++) + { + mchain[i].prev = seg[i].prev; + mchain[i].next = seg[i].next; + mchain[i].vnum = i; + vert[i].pt = seg[i].v0; + vert[i].vnext[0] = seg[i].next; /* next vertex */ + vert[i].vpos[0] = i; /* locn. of next vertex */ + vert[i].nextfree = 1; + } + + chain_idx = n; + mon_idx = 0; + mon[0] = 1; /* position of any vertex in the first */ + /* chain */ + +#endif + + /* traverse the polygon */ + if (tr[tr_start].u0 > 0) + traverse_polygon(0, tr_start, tr[tr_start].u0, TR_FROM_UP); + else if (tr[tr_start].d0 > 0) + traverse_polygon(0, tr_start, tr[tr_start].d0, TR_FROM_DN); + + /* return the number of polygons created */ + return newmon(); +} + + +/* recursively visit all the trapezoids */ +static int traverse_polygon(mcur, trnum, from, dir) + int mcur; + int trnum; + int from; + int dir; +{ + trap_t *t = &tr[trnum]; + int howsplit, mnew; + int v0, v1, v0next, v1next; + int retval, tmp; + int do_switch = FALSE; + + if ((trnum <= 0) || visited[trnum]) + return 0; + + visited[trnum] = TRUE; + + /* We have much more information available here. */ + /* rseg: goes upwards */ + /* lseg: goes downwards */ + + /* Initially assume that dir = TR_FROM_DN (from the left) */ + /* Switch v0 and v1 if necessary afterwards */ + + + /* special cases for triangles with cusps at the opposite ends. */ + /* take care of this first */ + if ((t->u0 <= 0) && (t->u1 <= 0)) + { + if ((t->d0 > 0) && (t->d1 > 0)) /* downward opening triangle */ + { + v0 = tr[t->d1].lseg; + v1 = t->lseg; + if (from == t->d1) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + } + } + else + { + retval = SP_NOSPLIT; /* Just traverse all neighbours */ + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + } + } + + else if ((t->d0 <= 0) && (t->d1 <= 0)) + { + if ((t->u0 > 0) && (t->u1 > 0)) /* upward opening triangle */ + { + v0 = t->rseg; + v1 = tr[t->u0].rseg; + if (from == t->u1) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + } + } + else + { + retval = SP_NOSPLIT; /* Just traverse all neighbours */ + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + } + } + + else if ((t->u0 > 0) && (t->u1 > 0)) + { + if ((t->d0 > 0) && (t->d1 > 0)) /* downward + upward cusps */ + { + v0 = tr[t->d1].lseg; + v1 = tr[t->u0].rseg; + retval = SP_2UP_2DN; + if (((dir == TR_FROM_DN) && (t->d1 == from)) || + ((dir == TR_FROM_UP) && (t->u1 == from))) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + } + } + else /* only downward cusp */ + { + if (_equal_to(&t->lo, &seg[t->lseg].v1)) + { + v0 = tr[t->u0].rseg; + v1 = seg[t->lseg].next; + + retval = SP_2UP_LEFT; + if ((dir == TR_FROM_UP) && (t->u0 == from)) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + } + } + else + { + v0 = t->rseg; + v1 = tr[t->u0].rseg; + retval = SP_2UP_RIGHT; + if ((dir == TR_FROM_UP) && (t->u1 == from)) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + } + } + } + } + else if ((t->u0 > 0) || (t->u1 > 0)) /* no downward cusp */ + { + if ((t->d0 > 0) && (t->d1 > 0)) /* only upward cusp */ + { + if (_equal_to(&t->hi, &seg[t->lseg].v0)) + { + v0 = tr[t->d1].lseg; + v1 = t->lseg; + retval = SP_2DN_LEFT; + if (!((dir == TR_FROM_DN) && (t->d0 == from))) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + } + } + else + { + v0 = tr[t->d1].lseg; + v1 = seg[t->rseg].next; + + retval = SP_2DN_RIGHT; + if ((dir == TR_FROM_DN) && (t->d1 == from)) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + } + } + } + else /* no cusp */ + { + if (_equal_to(&t->hi, &seg[t->lseg].v0) && + _equal_to(&t->lo, &seg[t->rseg].v0)) + { + v0 = t->rseg; + v1 = t->lseg; + retval = SP_SIMPLE_LRDN; + if (dir == TR_FROM_UP) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + } + } + else if (_equal_to(&t->hi, &seg[t->rseg].v1) && + _equal_to(&t->lo, &seg[t->lseg].v1)) + { + v0 = seg[t->rseg].next; + v1 = seg[t->lseg].next; + + retval = SP_SIMPLE_LRUP; + if (dir == TR_FROM_UP) + { + do_switch = TRUE; + mnew = make_new_monotone_poly(mcur, v1, v0); + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP); + } + else + { + mnew = make_new_monotone_poly(mcur, v0, v1); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN); + } + } + else /* no split possible */ + { + retval = SP_NOSPLIT; + traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP); + traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN); + traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP); + } + } + } + + return retval; +} + + +/* For each monotone polygon, find the ymax and ymin (to determine the */ +/* two y-monotone chains) and pass on this monotone polygon for greedy */ +/* triangulation. */ +/* Take care not to triangulate duplicate monotone polygons */ + +int triangulate_monotone_polygons(nvert, nmonpoly, op) + int nvert; + int nmonpoly; + int op[][3]; +{ + register int i; + point_t ymax, ymin; + int p, vfirst, posmax, posmin, v; + int vcount, processed; + +#ifdef DEBUG + for (i = 0; i < nmonpoly; i++) + { + fprintf(stderr, "\n\nPolygon %d: ", i); + vfirst = mchain[mon[i]].vnum; + p = mchain[mon[i]].next; + fprintf (stderr, "%d ", mchain[mon[i]].vnum); + while (mchain[p].vnum != vfirst) + { + fprintf(stderr, "%d ", mchain[p].vnum); + p = mchain[p].next; + } + } + fprintf(stderr, "\n"); +#endif + + op_idx = 0; + for (i = 0; i < nmonpoly; i++) + { + vcount = 1; + processed = FALSE; + vfirst = mchain[mon[i]].vnum; + ymax = ymin = vert[vfirst].pt; + posmax = posmin = mon[i]; + mchain[mon[i]].marked = TRUE; + p = mchain[mon[i]].next; + while ((v = mchain[p].vnum) != vfirst) + { + if (mchain[p].marked) + { + processed = TRUE; + break; /* break from while */ + } + else + mchain[p].marked = TRUE; + + if (_greater_than(&vert[v].pt, &ymax)) + { + ymax = vert[v].pt; + posmax = p; + } + if (_less_than(&vert[v].pt, &ymin)) + { + ymin = vert[v].pt; + posmin = p; + } + p = mchain[p].next; + vcount++; + } + + if (processed) /* Go to next polygon */ + continue; + + if (vcount == 3) /* already a triangle */ + { + op[op_idx][0] = mchain[p].vnum; + op[op_idx][1] = mchain[mchain[p].next].vnum; + op[op_idx][2] = mchain[mchain[p].prev].vnum; + op_idx++; + } + else /* triangulate the polygon */ + { + v = mchain[mchain[posmax].next].vnum; + if (_equal_to(&vert[v].pt, &ymin)) + { /* LHS is a single line */ + triangulate_single_polygon(nvert, posmax, TRI_LHS, op); + } + else + triangulate_single_polygon(nvert, posmax, TRI_RHS, op); + } + } + +#ifdef DEBUG + for (i = 0; i < op_idx; i++) + fprintf(stderr, "tri #%d: (%d, %d, %d)\n", i, op[i][0], op[i][1], + op[i][2]); +#endif + return op_idx; +} + + +/* A greedy corner-cutting algorithm to triangulate a y-monotone + * polygon in O(n) time. + * Joseph O-Rourke, Computational Geometry in C. + */ +static int triangulate_single_polygon(nvert, posmax, side, op) + int nvert; + int posmax; + int side; + int op[][3]; +{ + register int v; + int rc[SEGSIZE], ri = 0; /* reflex chain */ + int endv, tmp, vpos; + + if (side == TRI_RHS) /* RHS segment is a single segment */ + { + rc[0] = mchain[posmax].vnum; + tmp = mchain[posmax].next; + rc[1] = mchain[tmp].vnum; + ri = 1; + + vpos = mchain[tmp].next; + v = mchain[vpos].vnum; + + if ((endv = mchain[mchain[posmax].prev].vnum) == 0) + endv = nvert; + } + else /* LHS is a single segment */ + { + tmp = mchain[posmax].next; + rc[0] = mchain[tmp].vnum; + tmp = mchain[tmp].next; + rc[1] = mchain[tmp].vnum; + ri = 1; + + vpos = mchain[tmp].next; + v = mchain[vpos].vnum; + + endv = mchain[posmax].vnum; + } + + while ((v != endv) || (ri > 1)) + { + if (ri > 0) /* reflex chain is non-empty */ + { + if (CROSS(vert[v].pt, vert[rc[ri - 1]].pt, + vert[rc[ri]].pt) > 0) + { /* convex corner: cut if off */ + op[op_idx][0] = rc[ri - 1]; + op[op_idx][1] = rc[ri]; + op[op_idx][2] = v; + op_idx++; + ri--; + } + else /* non-convex */ + { /* add v to the chain */ + ri++; + rc[ri] = v; + vpos = mchain[vpos].next; + v = mchain[vpos].vnum; + } + } + else /* reflex-chain empty: add v to the */ + { /* reflex chain and advance it */ + rc[++ri] = v; + vpos = mchain[vpos].next; + v = mchain[vpos].vnum; + } + } /* end-while */ + + /* reached the bottom vertex. Add in the triangle formed */ + op[op_idx][0] = rc[ri - 1]; + op[op_idx][1] = rc[ri]; + op[op_idx][2] = v; + op_idx++; + ri--; + + return 0; +} + + Index: maptool/seidel/misc_seidel.c =================================================================== --- maptool/seidel/misc_seidel.c (revision 0) +++ maptool/seidel/misc_seidel.c (revision 0) @@ -0,0 +1,153 @@ +#include +#include +#include "triangulate.h" + +#ifdef __STDC__ +extern double log2(double); +#else +extern double log2(); +#endif + +static int choose_idx; +static int permute[SEGSIZE]; + + +/* Generate a random permutation of the segments 1..n */ +int generate_random_ordering(n) + int n; +{ + struct timeval tval; + struct timezone tzone; + register int i; + int m, st[SEGSIZE], *p; + + choose_idx = 1; + gettimeofday(&tval, &tzone); + srand48(tval.tv_sec); + + for (i = 0; i <= n; i++) + st[i] = i; + + p = st; + for (i = 1; i <= n; i++, p++) + { + m = lrand48() % (n + 1 - i) + 1; + permute[i] = p[m]; + if (m != 1) + p[m] = p[1]; + } + return 0; +} + + +/* Return the next segment in the generated random ordering of all the */ +/* segments in S */ +int choose_segment() +{ + int i; + +#ifdef DEBUG + fprintf(stderr, "choose_segment: %d\n", permute[choose_idx]); +#endif + return permute[choose_idx++]; +} + + +#ifdef STANDALONE + +/* Read in the list of vertices from infile */ +int read_segments(filename, genus) + char *filename; + int *genus; +{ + FILE *infile; + int ccount; + register int i; + int ncontours, npoints, first, last; + + if ((infile = fopen(filename, "r")) == NULL) + { + perror(filename); + return -1; + } + + fscanf(infile, "%d", &ncontours); + if (ncontours <= 0) + return -1; + + /* For every contour, read in all the points for the contour. The */ + /* outer-most contour is read in first (points specified in */ + /* anti-clockwise order). Next, the inner contours are input in */ + /* clockwise order */ + + ccount = 0; + i = 1; + + while (ccount < ncontours) + { + int j; + + fscanf(infile, "%d", &npoints); + first = i; + last = first + npoints - 1; + for (j = 0; j < npoints; j++, i++) + { + fscanf(infile, "%lf%lf", &seg[i].v0.x, &seg[i].v0.y); + if (i == last) + { + seg[i].next = first; + seg[i].prev = i-1; + seg[i-1].v1 = seg[i].v0; + } + else if (i == first) + { + seg[i].next = i+1; + seg[i].prev = last; + seg[last].v1 = seg[i].v0; + } + else + { + seg[i].prev = i-1; + seg[i].next = i+1; + seg[i-1].v1 = seg[i].v0; + } + + seg[i].is_inserted = FALSE; + } + + ccount++; + } + + *genus = ncontours - 1; + return i-1; +} + +#endif + + +/* Get log*n for given n */ +int math_logstar_n(n) + int n; +{ + register int i; + double v; + + for (i = 0, v = (double) n; v >= 1; i++) + v = log2(v); + + return (i - 1); +} + + +int math_N(n, h) + int n; + int h; +{ + register int i; + double v; + + for (i = 0, v = (int) n; i < h; i++) + v = log2(v); + + return (int) ceil((double) 1.0*n/v); +} Index: maptool/seidel/construct.c =================================================================== --- maptool/seidel/construct.c (revision 0) +++ maptool/seidel/construct.c (revision 0) @@ -0,0 +1,1060 @@ +#include "triangulate.h" +#include + + +node_t qs[QSIZE]; /* Query structure */ +trap_t tr[TRSIZE]; /* Trapezoid structure */ +segment_t seg[SEGSIZE]; /* Segment table */ + +static int q_idx; +static int tr_idx; + +/* Return a new node to be added into the query tree */ +static int newnode() +{ + if (q_idx < QSIZE) + return q_idx++; + else + { + fprintf(stderr, "newnode: Query-table overflow\n"); + return -1; + } +} + +/* Return a free trapezoid */ +static int newtrap() +{ + if (tr_idx < TRSIZE) + { + tr[tr_idx].lseg = -1; + tr[tr_idx].rseg = -1; + tr[tr_idx].state = ST_VALID; + return tr_idx++; + } + else + { + fprintf(stderr, "newtrap: Trapezoid-table overflow\n"); + return -1; + } +} + + +/* Return the maximum of the two points into the yval structure */ +static int _max(yval, v0, v1) + point_t *yval; + point_t *v0; + point_t *v1; +{ + if (v0->y > v1->y + C_EPS) + *yval = *v0; + else if (FP_EQUAL(v0->y, v1->y)) + { + if (v0->x > v1->x + C_EPS) + *yval = *v0; + else + *yval = *v1; + } + else + *yval = *v1; + + return 0; +} + + +/* Return the minimum of the two points into the yval structure */ +static int _min(yval, v0, v1) + point_t *yval; + point_t *v0; + point_t *v1; +{ + if (v0->y < v1->y - C_EPS) + *yval = *v0; + else if (FP_EQUAL(v0->y, v1->y)) + { + if (v0->x < v1->x) + *yval = *v0; + else + *yval = *v1; + } + else + *yval = *v1; + + return 0; +} + + +int _greater_than(v0, v1) + point_t *v0; + point_t *v1; +{ + if (v0->y > v1->y + C_EPS) + return TRUE; + else if (v0->y < v1->y - C_EPS) + return FALSE; + else + return (v0->x > v1->x); +} + + +int _equal_to(v0, v1) + point_t *v0; + point_t *v1; +{ + return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x)); +} + +int _greater_than_equal_to(v0, v1) + point_t *v0; + point_t *v1; +{ + if (v0->y > v1->y + C_EPS) + return TRUE; + else if (v0->y < v1->y - C_EPS) + return FALSE; + else + return (v0->x >= v1->x); +} + +int _less_than(v0, v1) + point_t *v0; + point_t *v1; +{ + if (v0->y < v1->y - C_EPS) + return TRUE; + else if (v0->y > v1->y + C_EPS) + return FALSE; + else + return (v0->x < v1->x); +} + + +/* Initilialise the query structure (Q) and the trapezoid table (T) + * when the first segment is added to start the trapezoidation. The + * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes + * + * 4 + * ----------------------------------- + * \ + * 1 \ 2 + * \ + * ----------------------------------- + * 3 + */ + +static int init_query_structure(segnum) + int segnum; +{ + int i1, i2, i3, i4, i5, i6, i7, root; + int t1, t2, t3, t4; + segment_t *s = &seg[segnum]; + + q_idx = tr_idx = 1; + memset((void *)tr, 0, sizeof(tr)); + memset((void *)qs, 0, sizeof(qs)); + + i1 = newnode(); + qs[i1].nodetype = T_Y; + _max(&qs[i1].yval, &s->v0, &s->v1); /* root */ + root = i1; + + qs[i1].right = i2 = newnode(); + qs[i2].nodetype = T_SINK; + qs[i2].parent = i1; + + qs[i1].left = i3 = newnode(); + qs[i3].nodetype = T_Y; + _min(&qs[i3].yval, &s->v0, &s->v1); /* root */ + qs[i3].parent = i1; + + qs[i3].left = i4 = newnode(); + qs[i4].nodetype = T_SINK; + qs[i4].parent = i3; + + qs[i3].right = i5 = newnode(); + qs[i5].nodetype = T_X; + qs[i5].segnum = segnum; + qs[i5].parent = i3; + + qs[i5].left = i6 = newnode(); + qs[i6].nodetype = T_SINK; + qs[i6].parent = i5; + + qs[i5].right = i7 = newnode(); + qs[i7].nodetype = T_SINK; + qs[i7].parent = i5; + + t1 = newtrap(); /* middle left */ + t2 = newtrap(); /* middle right */ + t3 = newtrap(); /* bottom-most */ + t4 = newtrap(); /* topmost */ + + tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval; + tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval; + tr[t4].hi.y = (double) (INFINITY); + tr[t4].hi.x = (double) (INFINITY); + tr[t3].lo.y = (double) -1* (INFINITY); + tr[t3].lo.x = (double) -1* (INFINITY); + tr[t1].rseg = tr[t2].lseg = segnum; + tr[t1].u0 = tr[t2].u0 = t4; + tr[t1].d0 = tr[t2].d0 = t3; + tr[t4].d0 = tr[t3].u0 = t1; + tr[t4].d1 = tr[t3].u1 = t2; + + tr[t1].sink = i6; + tr[t2].sink = i7; + tr[t3].sink = i4; + tr[t4].sink = i2; + + tr[t1].state = tr[t2].state = ST_VALID; + tr[t3].state = tr[t4].state = ST_VALID; + + qs[i2].trnum = t4; + qs[i4].trnum = t3; + qs[i6].trnum = t1; + qs[i7].trnum = t2; + + s->is_inserted = TRUE; + return root; +} + + +/* Retun TRUE if the vertex v is to the left of line segment no. + * segnum. Takes care of the degenerate cases when both the vertices + * have the same y--cood, etc. + */ + +static int is_left_of(segnum, v) + int segnum; + point_t *v; +{ + segment_t *s = &seg[segnum]; + double area; + + if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */ + { + if (FP_EQUAL(s->v1.y, v->y)) + { + if (v->x < s->v1.x) + area = 1.0; + else + area = -1.0; + } + else if (FP_EQUAL(s->v0.y, v->y)) + { + if (v->x < s->v0.x) + area = 1.0; + else + area = -1.0; + } + else + area = CROSS(s->v0, s->v1, (*v)); + } + else /* v0 > v1 */ + { + if (FP_EQUAL(s->v1.y, v->y)) + { + if (v->x < s->v1.x) + area = 1.0; + else + area = -1.0; + } + else if (FP_EQUAL(s->v0.y, v->y)) + { + if (v->x < s->v0.x) + area = 1.0; + else + area = -1.0; + } + else + area = CROSS(s->v1, s->v0, (*v)); + } + + if (area > 0.0) + return TRUE; + else + return FALSE; +} + + + +/* Returns true if the corresponding endpoint of the given segment is */ +/* already inserted into the segment tree. Use the simple test of */ +/* whether the segment which shares this endpoint is already inserted */ + +static int inserted(segnum, whichpt) + int segnum; + int whichpt; +{ + if (whichpt == FIRSTPT) + return seg[seg[segnum].prev].is_inserted; + else + return seg[seg[segnum].next].is_inserted; +} + +/* This is query routine which determines which trapezoid does the + * point v lie in. The return value is the trapezoid number. + */ + +int locate_endpoint(v, vo, r) + point_t *v; + point_t *vo; + int r; +{ + node_t *rptr = &qs[r]; + + switch (rptr->nodetype) + { + case T_SINK: + return rptr->trnum; + + case T_Y: + if (_greater_than(v, &rptr->yval)) /* above */ + return locate_endpoint(v, vo, rptr->right); + else if (_equal_to(v, &rptr->yval)) /* the point is already */ + { /* inserted. */ + if (_greater_than(vo, &rptr->yval)) /* above */ + return locate_endpoint(v, vo, rptr->right); + else + return locate_endpoint(v, vo, rptr->left); /* below */ + } + else + return locate_endpoint(v, vo, rptr->left); /* below */ + + case T_X: + if (_equal_to(v, &seg[rptr->segnum].v0) || + _equal_to(v, &seg[rptr->segnum].v1)) + { + if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */ + { + if (vo->x < v->x) + return locate_endpoint(v, vo, rptr->left); /* left */ + else + return locate_endpoint(v, vo, rptr->right); /* right */ + } + + else if (is_left_of(rptr->segnum, vo)) + return locate_endpoint(v, vo, rptr->left); /* left */ + else + return locate_endpoint(v, vo, rptr->right); /* right */ + } + else if (is_left_of(rptr->segnum, v)) + return locate_endpoint(v, vo, rptr->left); /* left */ + else + return locate_endpoint(v, vo, rptr->right); /* right */ + + default: + fprintf(stderr, "Haggu !!!!!\n"); + break; + } +} + + +/* Thread in the segment into the existing trapezoidation. The + * limiting trapezoids are given by tfirst and tlast (which are the + * trapezoids containing the two endpoints of the segment. Merges all + * possible trapezoids which flank this segment and have been recently + * divided because of its insertion + */ + +static int merge_trapezoids(segnum, tfirst, tlast, side) + int segnum; + int tfirst; + int tlast; + int side; +{ + int t, tnext, cond; + int ptnext; + + /* First merge polys on the LHS */ + t = tfirst; + while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo)) + { + if (side == S_LEFT) + cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) || + (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum))); + else + cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) || + (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum))); + + if (cond) + { + if ((tr[t].lseg == tr[tnext].lseg) && + (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */ + { /* merge them */ + /* Use the upper node as the new node i.e. t */ + + ptnext = qs[tr[tnext].sink].parent; + + if (qs[ptnext].left == tr[tnext].sink) + qs[ptnext].left = tr[t].sink; + else + qs[ptnext].right = tr[t].sink; /* redirect parent */ + + + /* Change the upper neighbours of the lower trapezoids */ + + if ((tr[t].d0 = tr[tnext].d0) > 0) + if (tr[tr[t].d0].u0 == tnext) + tr[tr[t].d0].u0 = t; + else if (tr[tr[t].d0].u1 == tnext) + tr[tr[t].d0].u1 = t; + + if ((tr[t].d1 = tr[tnext].d1) > 0) + if (tr[tr[t].d1].u0 == tnext) + tr[tr[t].d1].u0 = t; + else if (tr[tr[t].d1].u1 == tnext) + tr[tr[t].d1].u1 = t; + + tr[t].lo = tr[tnext].lo; + tr[tnext].state = ST_INVALID; /* invalidate the lower */ + /* trapezium */ + } + else /* not good neighbours */ + t = tnext; + } + else /* do not satisfy the outer if */ + t = tnext; + + } /* end-while */ + + return 0; +} + + +/* Add in the new segment into the trapezoidation and update Q and T + * structures. First locate the two endpoints of the segment in the + * Q-structure. Then start from the topmost trapezoid and go down to + * the lower trapezoid dividing all the trapezoids in between . + */ + +static int add_segment(segnum) + int segnum; +{ + segment_t s; + segment_t *so = &seg[segnum]; + int tu, tl, sk, tfirst, tlast, tnext; + int tfirstr, tlastr, tfirstl, tlastl; + int i1, i2, t, t1, t2, tn; + point_t tpt; + int tritop = 0, tribot = 0, is_swapped = 0; + int tmptriseg; + + s = seg[segnum]; + if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */ + { + int tmp; + tpt = s.v0; + s.v0 = s.v1; + s.v1 = tpt; + tmp = s.root0; + s.root0 = s.root1; + s.root1 = tmp; + is_swapped = TRUE; + } + + if ((is_swapped) ? !inserted(segnum, LASTPT) : + !inserted(segnum, FIRSTPT)) /* insert v0 in the tree */ + { + int tmp_d; + + tu = locate_endpoint(&s.v0, &s.v1, s.root0); + tl = newtrap(); /* tl is the new lower trapezoid */ + tr[tl].state = ST_VALID; + tr[tl] = tr[tu]; + tr[tu].lo.y = tr[tl].hi.y = s.v0.y; + tr[tu].lo.x = tr[tl].hi.x = s.v0.x; + tr[tu].d0 = tl; + tr[tu].d1 = 0; + tr[tl].u0 = tu; + tr[tl].u1 = 0; + + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + /* Now update the query structure and obtain the sinks for the */ + /* two trapezoids */ + + i1 = newnode(); /* Upper trapezoid sink */ + i2 = newnode(); /* Lower trapezoid sink */ + sk = tr[tu].sink; + + qs[sk].nodetype = T_Y; + qs[sk].yval = s.v0; + qs[sk].segnum = segnum; /* not really reqd ... maybe later */ + qs[sk].left = i2; + qs[sk].right = i1; + + qs[i1].nodetype = T_SINK; + qs[i1].trnum = tu; + qs[i1].parent = sk; + + qs[i2].nodetype = T_SINK; + qs[i2].trnum = tl; + qs[i2].parent = sk; + + tr[tu].sink = i1; + tr[tl].sink = i2; + tfirst = tl; + } + else /* v0 already present */ + { /* Get the topmost intersecting trapezoid */ + tfirst = locate_endpoint(&s.v0, &s.v1, s.root0); + tritop = 1; + } + + + if ((is_swapped) ? !inserted(segnum, FIRSTPT) : + !inserted(segnum, LASTPT)) /* insert v1 in the tree */ + { + int tmp_d; + + tu = locate_endpoint(&s.v1, &s.v0, s.root1); + + tl = newtrap(); /* tl is the new lower trapezoid */ + tr[tl].state = ST_VALID; + tr[tl] = tr[tu]; + tr[tu].lo.y = tr[tl].hi.y = s.v1.y; + tr[tu].lo.x = tr[tl].hi.x = s.v1.x; + tr[tu].d0 = tl; + tr[tu].d1 = 0; + tr[tl].u0 = tu; + tr[tl].u1 = 0; + + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu)) + tr[tmp_d].u0 = tl; + if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu)) + tr[tmp_d].u1 = tl; + + /* Now update the query structure and obtain the sinks for the */ + /* two trapezoids */ + + i1 = newnode(); /* Upper trapezoid sink */ + i2 = newnode(); /* Lower trapezoid sink */ + sk = tr[tu].sink; + + qs[sk].nodetype = T_Y; + qs[sk].yval = s.v1; + qs[sk].segnum = segnum; /* not really reqd ... maybe later */ + qs[sk].left = i2; + qs[sk].right = i1; + + qs[i1].nodetype = T_SINK; + qs[i1].trnum = tu; + qs[i1].parent = sk; + + qs[i2].nodetype = T_SINK; + qs[i2].trnum = tl; + qs[i2].parent = sk; + + tr[tu].sink = i1; + tr[tl].sink = i2; + tlast = tu; + } + else /* v1 already present */ + { /* Get the lowermost intersecting trapezoid */ + tlast = locate_endpoint(&s.v1, &s.v0, s.root1); + tribot = 1; + } + + /* Thread the segment into the query tree creating a new X-node */ + /* First, split all the trapezoids which are intersected by s into */ + /* two */ + + t = tfirst; /* topmost trapezoid */ + + while ((t > 0) && + _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo)) + /* traverse from top to bot */ + { + int t_sav, tn_sav; + sk = tr[t].sink; + i1 = newnode(); /* left trapezoid sink */ + i2 = newnode(); /* right trapezoid sink */ + + qs[sk].nodetype = T_X; + qs[sk].segnum = segnum; + qs[sk].left = i1; + qs[sk].right = i2; + + qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */ + qs[i1].trnum = t; + qs[i1].parent = sk; + + qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */ + qs[i2].trnum = tn = newtrap(); + tr[tn].state = ST_VALID; + qs[i2].parent = sk; + + if (t == tfirst) + tfirstr = tn; + if (_equal_to(&tr[t].lo, &tr[tlast].lo)) + tlastr = tn; + + tr[tn] = tr[t]; + tr[t].sink = i1; + tr[tn].sink = i2; + t_sav = t; + tn_sav = tn; + + /* error */ + + if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */ + { + fprintf(stderr, "add_segment: error\n"); + return 1; + //break; + } + + /* only one trapezoid below. partition t into two and make the */ + /* two resulting trapezoids t and tn as the upper neighbours of */ + /* the sole lower trapezoid */ + + else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0)) + { /* Only one trapezoid below */ + if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) + { /* continuation of a chain from abv. */ + if (tr[t].usave > 0) /* three upper neighbours */ + { + if (tr[t].uside == S_LEFT) + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = -1; + tr[tn].u1 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[tn].u0].d0 = tn; + tr[tr[tn].u1].d0 = tn; + } + else /* intersects in the right */ + { + tr[tn].u1 = -1; + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[t].u0; + tr[t].u0 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } + + tr[t].usave = tr[tn].usave = 0; + } + else /* No usave.... simple case */ + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d0 = tn; + } + } + else + { /* fresh seg. or upward cusp */ + int tmp_u = tr[t].u0; + int td0, td1; + if (((td0 = tr[tmp_u].d0) > 0) && + ((td1 = tr[tmp_u].d1) > 0)) + { /* upward cusp */ + if ((tr[td0].rseg > 0) && + !is_left_of(tr[td0].rseg, &s.v1)) + { + tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d1 = tn; + } + else /* cusp going leftwards */ + { + tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1; + tr[tr[t].u0].d0 = t; + } + } + else /* fresh segment */ + { + tr[tr[t].u0].d0 = t; + tr[tr[t].u0].d1 = tn; + } + } + + if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && + FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) + { /* bottom forms a triangle */ + + if (is_swapped) + tmptriseg = seg[segnum].prev; + else + tmptriseg = seg[segnum].next; + + if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0)) + { + /* L-R downward cusp */ + tr[tr[t].d0].u0 = t; + tr[tn].d0 = tr[tn].d1 = -1; + } + else + { + /* R-L downward cusp */ + tr[tr[tn].d0].u1 = tn; + tr[t].d0 = tr[t].d1 = -1; + } + } + else + { + if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0)) + { + if (tr[tr[t].d0].u0 == t) /* passes thru LHS */ + { + tr[tr[t].d0].usave = tr[tr[t].d0].u1; + tr[tr[t].d0].uside = S_LEFT; + } + else + { + tr[tr[t].d0].usave = tr[tr[t].d0].u0; + tr[tr[t].d0].uside = S_RIGHT; + } + } + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = tn; + } + + t = tr[t].d0; + } + + + else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0)) + { /* Only one trapezoid below */ + if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) + { /* continuation of a chain from abv. */ + if (tr[t].usave > 0) /* three upper neighbours */ + { + if (tr[t].uside == S_LEFT) + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = -1; + tr[tn].u1 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[tn].u0].d0 = tn; + tr[tr[tn].u1].d0 = tn; + } + else /* intersects in the right */ + { + tr[tn].u1 = -1; + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[t].u0; + tr[t].u0 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } + + tr[t].usave = tr[tn].usave = 0; + } + else /* No usave.... simple case */ + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d0 = tn; + } + } + else + { /* fresh seg. or upward cusp */ + int tmp_u = tr[t].u0; + int td0, td1; + if (((td0 = tr[tmp_u].d0) > 0) && + ((td1 = tr[tmp_u].d1) > 0)) + { /* upward cusp */ + if ((tr[td0].rseg > 0) && + !is_left_of(tr[td0].rseg, &s.v1)) + { + tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d1 = tn; + } + else + { + tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1; + tr[tr[t].u0].d0 = t; + } + } + else /* fresh segment */ + { + tr[tr[t].u0].d0 = t; + tr[tr[t].u0].d1 = tn; + } + } + + if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && + FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) + { /* bottom forms a triangle */ + int tmpseg; + + if (is_swapped) + tmptriseg = seg[segnum].prev; + else + tmptriseg = seg[segnum].next; + + if ((tmpseg > 0) && is_left_of(tmpseg, &s.v0)) + { + /* L-R downward cusp */ + tr[tr[t].d1].u0 = t; + tr[tn].d0 = tr[tn].d1 = -1; + } + else + { + /* R-L downward cusp */ + tr[tr[tn].d1].u1 = tn; + tr[t].d0 = tr[t].d1 = -1; + } + } + else + { + if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0)) + { + if (tr[tr[t].d1].u0 == t) /* passes thru LHS */ + { + tr[tr[t].d1].usave = tr[tr[t].d1].u1; + tr[tr[t].d1].uside = S_LEFT; + } + else + { + tr[tr[t].d1].usave = tr[tr[t].d1].u0; + tr[tr[t].d1].uside = S_RIGHT; + } + } + tr[tr[t].d1].u0 = t; + tr[tr[t].d1].u1 = tn; + } + + t = tr[t].d1; + } + + /* two trapezoids below. Find out which one is intersected by */ + /* this segment and proceed down that one */ + + else + { + int tmpseg = tr[tr[t].d0].rseg; + double y0, yt; + point_t tmppt; + int tnext, i_d0, i_d1; + + i_d0 = i_d1 = FALSE; + if (FP_EQUAL(tr[t].lo.y, s.v0.y)) + { + if (tr[t].lo.x > s.v0.x) + i_d0 = TRUE; + else + i_d1 = TRUE; + } + else + { + tmppt.y = y0 = tr[t].lo.y; + yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y); + tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x); + + if (_less_than(&tmppt, &tr[t].lo)) + i_d0 = TRUE; + else + i_d1 = TRUE; + } + + /* check continuity from the top so that the lower-neighbour */ + /* values are properly filled for the upper trapezoid */ + + if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) + { /* continuation of a chain from abv. */ + if (tr[t].usave > 0) /* three upper neighbours */ + { + if (tr[t].uside == S_LEFT) + { + tr[tn].u0 = tr[t].u1; + tr[t].u1 = -1; + tr[tn].u1 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[tn].u0].d0 = tn; + tr[tr[tn].u1].d0 = tn; + } + else /* intersects in the right */ + { + tr[tn].u1 = -1; + tr[tn].u0 = tr[t].u1; + tr[t].u1 = tr[t].u0; + tr[t].u0 = tr[t].usave; + + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } + + tr[t].usave = tr[tn].usave = 0; + } + else /* No usave.... simple case */ + { + tr[tn].u0 = tr[t].u1; + tr[tn].u1 = -1; + tr[t].u1 = -1; + tr[tr[tn].u0].d0 = tn; + } + } + else + { /* fresh seg. or upward cusp */ + int tmp_u = tr[t].u0; + int td0, td1; + if (((td0 = tr[tmp_u].d0) > 0) && + ((td1 = tr[tmp_u].d1) > 0)) + { /* upward cusp */ + if ((tr[td0].rseg > 0) && + !is_left_of(tr[td0].rseg, &s.v1)) + { + tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1; + tr[tr[tn].u0].d1 = tn; + } + else + { + tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1; + tr[tr[t].u0].d0 = t; + } + } + else /* fresh segment */ + { + tr[tr[t].u0].d0 = t; + tr[tr[t].u0].d1 = tn; + } + } + + if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && + FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) + { + /* this case arises only at the lowest trapezoid.. i.e. + tlast, if the lower endpoint of the segment is + already inserted in the structure */ + + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = -1; + tr[tr[t].d1].u0 = tn; + tr[tr[t].d1].u1 = -1; + + tr[tn].d0 = tr[t].d1; + tr[t].d1 = tr[tn].d1 = -1; + + tnext = tr[t].d1; + } + else if (i_d0) + /* intersecting d0 */ + { + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = tn; + tr[tr[t].d1].u0 = tn; + tr[tr[t].d1].u1 = -1; + + /* new code to determine the bottom neighbours of the */ + /* newly partitioned trapezoid */ + + tr[t].d1 = -1; + + tnext = tr[t].d0; + } + else /* intersecting d1 */ + { + tr[tr[t].d0].u0 = t; + tr[tr[t].d0].u1 = -1; + tr[tr[t].d1].u0 = t; + tr[tr[t].d1].u1 = tn; + + /* new code to determine the bottom neighbours of the */ + /* newly partitioned trapezoid */ + + tr[tn].d0 = tr[t].d1; + tr[tn].d1 = -1; + + tnext = tr[t].d1; + } + + t = tnext; + } + + tr[t_sav].rseg = tr[tn_sav].lseg = segnum; + } /* end-while */ + + /* Now combine those trapezoids which share common segments. We can */ + /* use the pointers to the parent to connect these together. This */ + /* works only because all these new trapezoids have been formed */ + /* due to splitting by the segment, and hence have only one parent */ + + tfirstl = tfirst; + tlastl = tlast; + merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT); + merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT); + + seg[segnum].is_inserted = TRUE; + return 0; +} + + +/* Update the roots stored for each of the endpoints of the segment. + * This is done to speed up the location-query for the endpoint when + * the segment is inserted into the trapezoidation subsequently + */ +static int find_new_roots(segnum) + int segnum; +{ + segment_t *s = &seg[segnum]; + + if (s->is_inserted) + return 0; + + s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0); + s->root0 = tr[s->root0].sink; + + s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1); + s->root1 = tr[s->root1].sink; + return 0; +} + + +/* Main routine to perform trapezoidation */ +int construct_trapezoids(nseg) + int nseg; +{ + register int i; + int root, h; + + /* Add the first segment and get the query structure and trapezoid */ + /* list initialised */ + + root = init_query_structure(choose_segment()); + + for (i = 1; i <= nseg; i++) + seg[i].root0 = seg[i].root1 = root; + + for (h = 1; h <= math_logstar_n(nseg); h++) + { + for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++) + if(add_segment(choose_segment())) { + return 1; + } + + /* Find a new root for each of the segment endpoints */ + for (i = 1; i <= nseg; i++) + find_new_roots(i); + } + + for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++) + if(add_segment(choose_segment())) { + return 1; + } + + return 0; +} + + Index: maptool/seidel/triangulate.h =================================================================== --- maptool/seidel/triangulate.h (revision 0) +++ maptool/seidel/triangulate.h (revision 0) @@ -0,0 +1,159 @@ +#ifndef _triangulate_h +#define _triangulate_h + +#include +#include +#include +#include "types.h" + +typedef struct { + double x, y; +} point_t, vector_t; + + +/* Segment attributes */ + +typedef struct { + point_t v0, v1; /* two endpoints */ + int is_inserted; /* inserted in trapezoidation yet ? */ + int root0, root1; /* root nodes in Q */ + int next; /* Next logical segment */ + int prev; /* Previous segment */ +} segment_t; + + +/* Trapezoid attributes */ + +typedef struct { + int lseg, rseg; /* two adjoining segments */ + point_t hi, lo; /* max/min y-values */ + int u0, u1; + int d0, d1; + int sink; /* pointer to corresponding in Q */ + int usave, uside; /* I forgot what this means */ + int state; +} trap_t; + + +/* Node attributes for every node in the query structure */ + +typedef struct { + int nodetype; /* Y-node or S-node */ + int segnum; + point_t yval; + int trnum; + int parent; /* doubly linked DAG */ + int left, right; /* children */ +} node_t; + + +typedef struct { + int vnum; + int next; /* Circularly linked list */ + int prev; /* describing the monotone */ + int marked; /* polygon */ +} monchain_t; + + +typedef struct { + point_t pt; + int vnext[4]; /* next vertices for the 4 chains */ + int vpos[4]; /* position of v in the 4 chains */ + int nextfree; +} vertexchain_t; + + +/* Node types */ + +#define T_X 1 +#define T_Y 2 +#define T_SINK 3 + + +#define SEGSIZE 1000000 /* max# of segments. Determines how */ + /* many points can be specified as */ + /* input. If your datasets have large */ + /* number of points, increase this */ + /* value accordingly. */ + +#define QSIZE 8*SEGSIZE /* maximum table sizes */ +#define TRSIZE 4*SEGSIZE /* max# trapezoids */ + +#define TRUE 1 +#define FALSE 0 + +#define FIRSTPT 1 /* checking whether pt. is inserted */ +#define LASTPT 2 + + +#define INFINITY 1<<30 +//#define C_EPS 1.0e-7 /* tolerance value: Used for making */ +#define C_EPS 1.0e-7 /* tolerance value: Used for making */ + /* all decisions about collinearity or */ + /* left/right of segment. Decrease */ + /* this value if the input points are */ + /* spaced very close together */ + + +#define S_LEFT 1 /* for merge-direction */ +#define S_RIGHT 2 + + +#define ST_VALID 1 /* for trapezium state */ +#define ST_INVALID 2 + + +#define SP_SIMPLE_LRUP 1 /* for splitting trapezoids */ +#define SP_SIMPLE_LRDN 2 +#define SP_2UP_2DN 3 +#define SP_2UP_LEFT 4 +#define SP_2UP_RIGHT 5 +#define SP_2DN_LEFT 6 +#define SP_2DN_RIGHT 7 +#define SP_NOSPLIT -1 + +#define TR_FROM_UP 1 /* for traverse-direction */ +#define TR_FROM_DN 2 + +#define TRI_LHS 1 +#define TRI_RHS 2 + + +#define MAX(a, b) (((a) > (b)) ? (a) : (b)) +#define MIN(a, b) (((a) < (b)) ? (a) : (b)) + +#define CROSS(v0, v1, v2) (((v1).x - (v0).x)*((v2).y - (v0).y) - \ + ((v1).y - (v0).y)*((v2).x - (v0).x)) + +#define DOT(v0, v1) ((v0).x * (v1).x + (v0).y * (v1).y) + +#define FP_EQUAL(s, t) (fabs(s - t) <= C_EPS) + + + +/* Global variables */ + +extern node_t qs[QSIZE]; /* Query structure */ +extern trap_t tr[TRSIZE]; /* Trapezoid structure */ +extern segment_t seg[SEGSIZE]; /* Segment table */ + + +/* Functions */ + +extern int monotonate_trapezoids(int); +extern int triangulate_monotone_polygons(int, int, int (*)[3]); + +extern int _greater_than(point_t *, point_t *); +extern int _equal_to(point_t *, point_t *); +extern int _greater_than_equal_to(point_t *, point_t *); +extern int _less_than(point_t *, point_t *); +extern int locate_endpoint(point_t *, point_t *, int); +extern int construct_trapezoids(int); + +extern int generate_random_ordering(int); +extern int choose_segment(void); +extern int read_segments(char *, int *); +extern int math_logstar_n(int); +extern int math_N(int, int); + +#endif /* triangulate_h */ Index: maptool/seidel/interface.h =================================================================== --- maptool/seidel/interface.h (revision 0) +++ maptool/seidel/interface.h (revision 0) @@ -0,0 +1,12 @@ +#ifndef __interface_h +#define __interface_h + +#include "types.h" + +#define TRUE 1 +#define FALSE 0 + +extern int triangulate_polygon(s32, s32*, s32 (*)[2], s32 (*)[3]); +extern int is_point_inside_polygon(double *); + +#endif /* __interface_h */ Index: maptool/geom.c =================================================================== --- maptool/geom.c (revision 5005) +++ maptool/geom.c (working copy) @@ -19,6 +19,7 @@ #include #include #include "maptool.h" +#include "transform.h" void geom_coord_copy(struct coord *from, struct coord *to, int count, int reverse) @@ -357,6 +358,127 @@ return ret; } +/* + * create groups of segments by merging segments with common endpoint + */ +GList* geom_poly_segments_group(GList *in, GList *out) +{ + out = g_list_copy(in); + + int changed; + do { + changed = 0; + GList*l1 = out; + while(l1) { + struct geom_poly_segment *l1_ = l1->data; + int size_1 = (l1_->last-l1_->first)+1; + if(((l1_->type!=geom_poly_segment_type_way_outer) && (l1_->type!=geom_poly_segment_type_way_inner)) || (size_1<1)) { + GList*l1_d = l1; + l1 = g_list_next(l1); + if(!l1) { + break; + } + l1_ = l1->data; + out = g_list_remove(out, l1_d->data); + continue; + } + GList*l2 = l1; + while(l2) { + struct geom_poly_segment *l2_ = l2->data; + int size_2 = (l2_->last-l2_->first)+1; + if(((l2_->type!=geom_poly_segment_type_way_outer) && (l2_->type!=geom_poly_segment_type_way_inner)) || (size_2<1)) { + GList*l2_d = l2; + l2 = g_list_next(l2); + if(!l2) { + break; + } + l2_ = l2->data; + out = g_list_remove(out, l2_d->data); + continue; + } + //suppose that inner segments are only connected with other inner segments + if( l1!=l2 && l1_ && l2_ && l1_->first && l2_->last && + (coord_equal(l1_->first, l2_->first) || /* have_common_endpoints(l1,l2)*/ + coord_equal(l1_->first, l2_->last) || + coord_equal(l1_->last, l2_->first) || + coord_equal(l1_->last, l2_->last)) && + !coord_equal(l1_->first,l1_->last) && !coord_equal(l2_->first,l2_->last) /*do not group closed polys*/ ) { + int reverse_1 = 0; + int reverse_2 = 0; + int order_12 = 1; + //merge coords of (l1,l2) (append coords to the front or rear, perhaps reversed) + //l1.front==l2.front + //l2 reversed ... l1 inorder + if (coord_equal(l1_->first,l2_->first)) { + reverse_1 = 0; + reverse_2 = 1; + order_12 = 0; + } + //l1.rear==l2.front + //l1 inorder ... l2 inorder + else if (coord_equal(l1_->last,l2_->first)) { + reverse_1 = 0; + reverse_2 = 0; + order_12 = 1; + } + //l1.front==l2.rear + //l2 inorder ... l1 inorder + else if (coord_equal(l1_->first,l2_->last)) { + reverse_1 = 0; + reverse_2 = 0; + order_12 = 0; + } + //l1.rear==l2.rear + //l1 inorder ... l2 reversed + else if (coord_equal(l1_->last,l2_->last)) { + reverse_1 = 0; + reverse_2 = 1; + order_12 = 1; + } + //reverse l1 or l2 if necessary + if (reverse_1) { + int ii; + for (ii=0;iifirst[ii]; + l1_->first[ii] = l1_->first[size_1-ii-1]; + l1_->first[size_1-ii-1] = tmp; + } + } + if (reverse_2) { + int ii; + for (ii=0;iifirst[ii]; + l2_->first[ii] = l2_->first[size_2-ii-1]; + l2_->first[size_2-ii-1] = tmp; + } + } + //append l1 to l2 or l2 to l1 appropriately (with the common point added only once) + if(order_12 && size_2>1) { + l1_->first = g_realloc(l1_->first,sizeof(struct coord)*(size_1+size_2)-1); + memcpy(l1_->first+size_1,l2_->first,(size_2-1)*sizeof(struct coord)); + l1_->last += size_2-1; + } else if( !order_12 && size_1>1) { + l2_->first = g_realloc(l2_->first,sizeof(struct coord)*(size_1+size_2-1)); + memcpy(l2_->first+size_2,l1_->first,(size_1-1)*sizeof(struct coord)); + l2_->last += size_1-1; + } + //remove l1 or l2 appropriately + if (order_12) { + out = g_list_remove(out, l2->data); + } else { + out = g_list_remove(out, l1->data); + } + changed=1; + } + l2 = g_list_next(l2); + } + l1 = g_list_next(l1); + } + } while(changed); +return out; +} + + int geom_poly_segments_point_inside(GList *in, struct coord *c) { @@ -620,3 +742,381 @@ item_bin_write_clipped(ib_in, param, out); } } + + + + +static int coord_dot_prod (struct coord* c1, struct coord* c2) +{ + return c1->x*c2->x + c1->y*c2->y; +} + + +static void coord_sub (struct coord* c1, struct coord* c2, struct coord* cout) +{ + cout->x = c1->x - c2->x; + cout->y = c1->y - c2->y; +} + +double DistanceSquared(struct coord* v0, struct coord* v1, struct coord*p) +{ + +/*handle the case of vertical line segments*/ + if (v0->x==v1->x && p->x==v0->x && + ( (v0->y<=p->y && p->y<=v1->y) || (v1->y<=p->y && p->y<=v0->y) ) + ) { + return 0.0; + } + + if (v0->y==v1->y && p->y==v0->y && + ( (v0->x<=p->x && p->x<=v1->x) || (v1->x<=p->x && p->x<=v0->x) ) + ) { + return 0.0; + } + + + double vx = v0->x - p->x, vy = v0->y - p->y, ux = v1->x - v0->x, uy = v1->y - v0->y; + double length = ux * ux + uy * uy; + + double det = (-vx * ux) + (-vy * uy); //if this is < 0 or > length then its outside the line segment + if (det < 0) { + return (v0->x - p->x) * (v0->x - p->x) + (v0->y - p->y) * (v0->y - p->y); + } + if (det > length) { + return (v1->x - p->x) * (v1->x - p->x) + (v1->y - p->y) * (v1->y - p->y); + } + + det = ux * vy - uy * vx; + return (det * det) / length; +} + + +//TODO merge intersect tests +int +self_intersect_test_segment(struct geom_poly_segment *gs) +{ + int vertices_count = (gs->last-gs->first)+1; + if(vertices_count<4) { + return 0; + } + struct coord*vertices = (struct coord*)gs->first; + if(coord_equal(&vertices[0],&vertices[vertices_count-1])) { + vertices_count--; + } + int i; + for (i = 0; i < vertices_count; ++i) { + if (i < vertices_count - 1) { + int h; + for (h = i + 1; h < vertices_count; ++h) + { + // Do two vertices lie on top of one another? + if ( i!=0 && h-i!=1 && h-i!=-1 && vertices[i].x == vertices[h].x && vertices[i].y == vertices[h].y ) { + return 1; + } + } + } + + int j = (i + 1) % vertices_count; + struct coord iToj; + coord_sub(&vertices[j], &vertices[i], &iToj); + struct coord iTojNormal; + iTojNormal.x = iToj.y; + iTojNormal.y = -iToj.x; + // i is the first vertex and j is the second + int startK = (j + 1) % vertices_count; + int endK = (i - 1 + vertices_count) % vertices_count; + endK += startK < endK ? 0 : startK + 1; + int k = startK; + struct coord iTok; + coord_sub(&vertices[k], &vertices[i], &iTok); + int onLeftSide = coord_dot_prod(&iTok, &iTojNormal) >= 0; + struct coord prevK = vertices[k]; + ++k; + for (; k <= endK; ++k) + { + int modK = k % vertices_count; + coord_sub(&vertices[modK], &vertices[i], &iTok); + if (onLeftSide != coord_dot_prod(&iTok, &iTojNormal) >= 0) + { + struct coord prevKtoK, idiff, jdiff; + coord_sub(&vertices[modK], &prevK, &prevKtoK); + struct coord prevKtoKNormal; + prevKtoKNormal.x = prevKtoK.y; + prevKtoKNormal.y = -prevKtoK.x; + coord_sub(&vertices[i],&prevK,&idiff); + coord_sub(&vertices[j],&prevK,&jdiff); + if ((coord_dot_prod(&idiff, &prevKtoKNormal) >= 0) != (coord_dot_prod(&jdiff, &prevKtoKNormal) >= 0)) + { + return 1; + } + } + onLeftSide = coord_dot_prod(&iTok, &iTojNormal) > 0; + prevK = vertices[modK]; + } + } + //check if a vertex lies on a lineseg defined by other vertices + struct coord* c = gs->first; + struct coord* c2; + while (c+1<=gs->last) { + c2 = gs->first; + while (c2<=gs->last) { + if ( !coord_equal(c,c2) && !coord_equal(c+1,c2) && DistanceSquared(c,c+1,c2)<.01) { + return 1; + } + c2++; + } + c++; + } + + if (!coord_equal(gs->first,gs->last)) { + c2 = gs->first; + while (c2<=gs->last) { + if ( !coord_equal(gs->first,c2) && !coord_equal(gs->last,c2) && DistanceSquared(gs->first,gs->last,c2)<.01) { + return 1; + } + c2++; + } + } + return 0; + +} + + +int self_intersect_test(struct item_bin*ib) +{ + int vertices_count = ib->clen/2; + if(vertices_count<4) { + return 0; + } + struct coord*vertices = (struct coord*)(ib+1); + + if(coord_equal(&vertices[0],&vertices[vertices_count-1])) { + vertices_count--; + } + int i; + for (i = 0; i < vertices_count; ++i) { + if (i < vertices_count - 1) { + int h; + for (h = i + 1; h < vertices_count; ++h) + { + // Do two vertices lie on top of one another? + if ( i!=0 && h-i!=1 && h-i!=-1 && vertices[i].x == vertices[h].x && vertices[i].y == vertices[h].y ) { + return 1; + } + } + } + + int j = (i + 1) % vertices_count; + struct coord iToj; + coord_sub(&vertices[j], &vertices[i], &iToj); + struct coord iTojNormal; + iTojNormal.x = iToj.y; + iTojNormal.y = -iToj.x; + // i is the first vertex and j is the second + int startK = (j + 1) % vertices_count; + int endK = (i - 1 + vertices_count) % vertices_count; + endK += startK < endK ? 0 : startK + 1; + int k = startK; + struct coord iTok; + coord_sub(&vertices[k], &vertices[i], &iTok); + int onLeftSide = coord_dot_prod(&iTok, &iTojNormal) >= 0; + struct coord prevK = vertices[k]; + ++k; + for (; k <= endK; ++k) + { + int modK = k % vertices_count; + coord_sub(&vertices[modK], &vertices[i], &iTok); + if (onLeftSide != coord_dot_prod(&iTok, &iTojNormal) >= 0) + { + struct coord prevKtoK, idiff, jdiff; + coord_sub(&vertices[modK], &prevK, &prevKtoK); + struct coord prevKtoKNormal; + prevKtoKNormal.x = prevKtoK.y; + prevKtoKNormal.y = -prevKtoK.x; + coord_sub(&vertices[i],&prevK,&idiff); + coord_sub(&vertices[j],&prevK,&jdiff); + if ((coord_dot_prod(&idiff, &prevKtoKNormal) >= 0) != (coord_dot_prod(&jdiff, &prevKtoKNormal) >= 0)) + { + if(1/*i!=j-1 && i!=j+1*/) { + return 1; + } + } + } + onLeftSide = coord_dot_prod(&iTok, &iTojNormal) > 0; + prevK = vertices[modK]; + } + } + return 0; +} + +/* + * returns true if coord c is inside poly segment poly + */ +int geom_point_in_polygon(struct coord*c, struct geom_poly_segment*poly) +{ + //TODO handle closed/opened polygon cases + int counter = 0; + int i,N; + double xinters; + struct coord p1,p2; + + if(!c || !poly) { + return 0; + } + N = poly->last-poly->first+1; + p1 = poly->first[0]; + + for (i=1;i<=N;i++) { + p2 = poly->first[i % N]; + if (c->y > MIN(p1.y,p2.y)) { + if (c->y <= MAX(p1.y,p2.y)) { + if (c->x <= MAX(p1.x,p2.x)) { + if (p1.y != p2.y) { + xinters = (c->y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x; + if (p1.x == p2.x || c->x <= xinters) + counter++; + } + } + } + } + p1 = p2; + } + + if (counter % 2 == 0) + return 0; + else + return 1; +} + +/* + * returns 1 if poly1 is approximately inside poly2 (every point of p1 is inside poly2) + */ +int geom_polygon_in_polygon(struct geom_poly_segment*poly1, struct geom_poly_segment*poly2) +{ + int ret=1; + struct coord *c = poly1->first; + do { + if(!geom_point_in_polygon(c,poly2)) { + ret=0; + break; + } + } while(c++ != poly1->last); + return ret; +} + + +int geom_polygon_is_clockwise(struct geom_poly_segment*poly) +{ + double sum=0; //sum of signed areas + struct coord *c = poly->first; + struct coord *c2; + do { + c2 = c+1; + sum += (double)c->x*c2->y-(double)c2->x*c->y; + } while(++c <= (poly->last-1)); + if(!coord_equal(poly->first,poly->last)) { + c =poly->last; + c2=poly->first; + sum += (double)c->x*c2->y-(double)c2->x*c->y; + } + return sum<0 ? 1 : 0; +} + +void geom_polygon_reverse(struct geom_poly_segment*poly) +{ + int n = poly->last-poly->first+1; + int ii; + for(ii=0;iifirst[ii]; + poly->first[ii] = poly->first[n-1-ii]; + poly->first[n-1-ii] = tmp; + } +} + +/* + * removes subsequent equivalent coords from segment + */ +void geom_polygon_unique(struct geom_poly_segment*poly) +{ + struct coord *c = poly->first; + struct coord *dst = NULL; + int unique_cnt=2; + int rearrange=0; + if(poly->last<=poly->first) { + return; + } + do { + struct coord *c2 = c+1; + if(!coord_equal(c,c2)) { + ++unique_cnt; + } else { + rearrange=1; + } + } while(++c != poly->last); + + if(rearrange) { + c = poly->first; + struct coord*new_segment_data = g_malloc(sizeof(struct coord)*unique_cnt); + dst=new_segment_data; + do { + struct coord *c2 = c+1; + *dst++ = *c; + while(c2<=poly->last && coord_equal(c,c2)) { + ++c; + ++c2; + } + } while(++c <= poly->last); + g_free(poly->first); + poly->first=new_segment_data; + poly->last=dst-1; + } +} + +int geom_triangle_ccw(struct coord*a,struct coord*b,struct coord*c) +{ + return (c->y-a->y)*(b->x-a->x) >= (b->y-a->y)*(c->x-a->x); +} + +int geom_lineseg_intersect(struct coord*a,struct coord*b,struct coord*c,struct coord*d) +{ + if(DistanceSquared(a,b,c)<.1 || DistanceSquared(a,b,d)<.1 || + DistanceSquared(c,d,a)<.1 || DistanceSquared(c,d,b)<.1 ) { + return 1; + } + + int d1 = transform_distance_line_sq(a,b,c,NULL); + int d2 = transform_distance_line_sq(a,b,d,NULL); + int d3 = transform_distance_line_sq(c,d,a,NULL); + int d4 = transform_distance_line_sq(c,d,b,NULL); + if(!d1 || !d2 || !d3 || !d4) { + return 1; + } + 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)); +} + +int geom_poly_intersect(struct geom_poly_segment*a,struct geom_poly_segment*b) +{ + struct coord*curr_a = a->first; + struct coord*curr_b = b->first; + while (curr_a+1 <= a->last) { //poly a segments + if(b->last!=b->first && geom_lineseg_intersect(curr_a,curr_a+1,b->last,b->first)) { + return 1; + break; + } + struct coord*curr_b = b->first; + while (curr_b+1 <= b->last) { //poly b segments + if(geom_lineseg_intersect(curr_a,curr_a+1,curr_b,curr_b+1)) { + return 1; + break; + } + curr_b++; + } + curr_a++; + } + if(b->last!=b->first && a->last!=a->first && geom_lineseg_intersect(a->last,a->first,b->last,b->first)) { + return 1; + } + + return 0; +} + Index: maptool/Makefile.am =================================================================== --- maptool/Makefile.am (revision 5005) +++ maptool/Makefile.am (working copy) @@ -5,6 +5,6 @@ endif AM_CPPFLAGS = @NAVIT_CFLAGS@ -I$(top_srcdir)/navit @ZLIB_CFLAGS@ @POSTGRESQL_CFLAGS@ -DMODULE=maptool -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 +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 seidel/interface.h seidel/triangulate.h seidel/misc_seidel.c seidel/construct.c seidel/monotone.c seidel/tri.c maptool_SOURCES = maptool.c maptool_LDADD = libmaptool.la ../libnavit.la @NAVIT_LIBS@ @WORDEXP_LIBS@ @ZLIB_LIBS@ @POSTGRESQL_LIBS@ @CRYPTO_LIBS@ @INTLLIBS@ @LIBC_LIBS@ Index: maptool/boundaries.c =================================================================== --- maptool/boundaries.c (revision 5005) +++ maptool/boundaries.c (working copy) @@ -19,6 +19,8 @@ #include #include #include "maptool.h" +#include "types.h" +#include "seidel/triangulate.h" char * osm_tag_value(struct item_bin *ib, char *key) @@ -70,6 +72,17 @@ while ((ib=read_item(boundaries))) { char *member=NULL; struct boundary *boundary=g_new0(struct boundary, 1); + boundary->area_item_types = NULL; + enum item_type* p_item_type; + if((p_item_type = item_bin_get_attr(ib,attr_item_types,NULL))) { + int i=0; + do { + boundary->area_item_types = g_realloc(boundary->area_item_types, (i+1)*sizeof(enum item_type)); + boundary->area_item_types[i] = *p_item_type; + i++; + boundary->area_item_types[i] = type_none; + } while(*p_item_type++!=type_none); + } char *admin_level=osm_tag_value(ib, "admin_level"); char *iso=osm_tag_value(ib, "ISO3166-1"); /* disable spain for now since it creates a too large index */ @@ -192,26 +205,266 @@ return g_list_prepend(list, boundary); } +//creates inner attr for the first outer +struct attr* +create_triangle_attr(GList*grouped_segments) +{ + struct attr *ret = NULL; + struct attr *triangle_attr = NULL; + GList *l = grouped_segments; + while (l) { + struct geom_poly_segment *gs=l->data; + if (coord_is_equal(*gs->first,*gs->last) && (gs->last-gs->first)>0) { + gs->last--; + } + if(gs->type==geom_poly_segment_type_way_outer && (gs->last-gs->first+1)>2 ) { //use the first outer + s32* polylen_buf = NULL; + s32 polylen = 0; + s32* vertices_buf = NULL; + GList *k = grouped_segments; + s32 prev_size = 0; + s32 coordsize_in_int=(sizeof(struct coord))/(sizeof(s32)); + if(geom_polygon_is_clockwise(gs)) { // outer polygons should be CCW + geom_polygon_reverse(gs); + } +#if 0 +printf("CRD:START\n"); +fflush(stdout); +#endif + while(k) { //loop through inner polygons + struct geom_poly_segment *gs2=k->data; + //if actual inner poly inside actual outer, add inner to the list of to be added inner polygons + if (coord_is_equal(*gs2->first,*gs2->last) && (gs2->last-gs2->first)>0) { + gs2->last--; + } + + 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)) { + //make sure that outer is CCW, inner is CW + if(!geom_polygon_is_clockwise(gs2)) { // inner polygons should be CW + geom_polygon_reverse(gs2); + } + if(!polylen_buf) { //collect info about the outer + polylen_buf = g_realloc(polylen_buf,sizeof(s32)*coordsize_in_int); + polylen_buf[0] = gs->last-gs->first+1; + ++polylen; + +#if 0 +struct coord*c=gs->first; +do { + printf("CRD:%d %d\n",c->x,c->y); + c++; +} while(c<=gs->last); +printf("CRD:\n"); +fflush(stdout); +#endif + vertices_buf = g_realloc(vertices_buf,sizeof(s32)*(coordsize_in_int*(gs->last-gs->first+1))); + prev_size += coordsize_in_int*(gs->last-gs->first+1); + memcpy(vertices_buf,gs->first,sizeof(s32)*coordsize_in_int*(gs->last-gs->first+1)); + } +#if 0 +struct coord*c=gs2->first; +do { + printf("CRD:%d %d\n",c->x,c->y); + c++; +} while(c<=gs2->last); +printf("CRD:\n"); +fflush(stdout); +#endif + + /* + */ + polylen_buf = g_realloc(polylen_buf,sizeof(s32)*coordsize_in_int*(++polylen)); + polylen_buf[polylen-1] = gs2->last-gs2->first+1; + vertices_buf = g_realloc(vertices_buf,sizeof(s32)*(prev_size+coordsize_in_int*(gs2->last-gs2->first+1))); + memcpy(&vertices_buf[prev_size],gs2->first,sizeof(s32)*coordsize_in_int*(gs2->last-gs2->first+1)); + prev_size+=coordsize_in_int*(gs2->last-gs2->first+1); + } + k=g_list_next(k); + } +#if 0 +printf("CRD:END\n"); +fflush(stdout); +#endif + +struct coord*c = (struct coord*)vertices_buf; +struct coord*startcoord = (struct coord*)vertices_buf; +int duplicate=0; +while(c<(startcoord+prev_size/coordsize_in_int) && duplicate==0) { + struct coord*c2 = c+1; + while(c2<(startcoord+prev_size/coordsize_in_int)) { + if(coord_is_equal(*c,*c2)) { + duplicate=1; + break; + } + c2++; + } + c++; +} +if( duplicate==0 && prev_size0 && polylen>0) { + int ii=0; + s32 (*triangles)[3] = g_malloc0(sizeof(s32)*tri_num*3); +#if 0 +printf("polylen %d\n",polylen); +for(ii=0;iitype = attr_triangles; + triangle_attr->u.data = pit; + + g_free(triangles); + } + + ret = triangle_attr; + + g_free(polylen_buf); + g_free(vertices_buf); + return ret; + } + g_free(polylen_buf); + g_free(vertices_buf); + } + +//TODO merge subsequent repeated points + + l=g_list_next(l); + } + return NULL; +} + +/* + * removes inner polygons intersecting the outer and other inners + * expects a grouped polygon lists + */ +GList* remove_intersecting_segments(GList*segments) +{ + //find outer + GList*curr=segments; + GList*curr2=segments; + struct geom_poly_segment*outer=NULL; + int changed; + int changed_g; + while(curr) { //find outer + struct geom_poly_segment *gs = curr->data; + changed=0; + if ( gs->last-gs->first<=0 || (gs->type==geom_poly_segment_type_way_outer && self_intersect_test_segment(gs)) ) { + curr=g_list_next(curr); + segments=g_list_remove(segments,gs); + changed=1; + } else if(gs->type==geom_poly_segment_type_way_outer) { + outer=gs; + } + if(!changed) { + curr = g_list_next(curr); + } + } + + curr=segments; + while(outer && curr) { + changed=0; + struct geom_poly_segment *gs = curr->data; + if(gs->type==geom_poly_segment_type_way_inner) { + if(geom_poly_intersect(outer,gs) || self_intersect_test_segment(gs) || geom_polygon_in_polygon(outer,gs) ) { + curr=g_list_next(curr); + segments=g_list_remove(segments,gs); + changed=1; + } + } else if(gs->type==geom_poly_segment_type_way_outer && gs!=outer) { + } + if(!changed) { + curr = g_list_next(curr); + } + } + + do { + changed_g=0; + curr=segments; + + while(curr) { + struct geom_poly_segment *gs = curr->data; + curr2 = g_list_next(curr); + while(curr2) { + struct geom_poly_segment *gs2 = curr2->data; + changed=0; + if(self_intersect_test_segment(gs2)) { + curr2=g_list_next(curr2); + segments=g_list_remove(segments,gs2); + changed=1; + changed_g=1; + continue; + } + if(gs->type==geom_poly_segment_type_way_inner && gs2->type==geom_poly_segment_type_way_inner) { + if( geom_polygon_in_polygon(gs2,gs) || /*currently we do not support nested inners*/ + geom_poly_intersect(gs,gs2)) { + + curr2=g_list_next(curr2); + segments=g_list_remove(segments,gs2); + changed=1; + changed_g=1; + } + } + if(!changed) { + curr2 = g_list_next(curr2); + } + } + curr = g_list_next(curr); + } + + } while(changed_g); + return segments; +} + static GList * process_boundaries_finish(GList *boundaries_list) { - GList *l,*sl,*l2,*ln; + GList *l,*sl; GList *ret=NULL; l=boundaries_list; while (l) { struct boundary *boundary=l->data; int first=1; - FILE *f=NULL,*fu=NULL; + FILE *f=NULL,*fu=NULL,*ways_split=NULL; + int to_write = 0; + int have_inner = 0; if (boundary->country) { char *name=g_strdup_printf("country_%s_poly",boundary->iso2); f=tempfile("",name,1); g_free(name); } + if (boundary->area_item_types) { + ways_split=tempfile("","ways_split",2); + struct item_bin *ib=item_bin_relation_area; + item_bin_init(ib, boundary->area_item_types[0]); + } boundary->sorted_segments=geom_poly_segments_sort(boundary->segments, geom_poly_segment_type_way_right_side); sl=boundary->sorted_segments; while (sl) { struct geom_poly_segment *gs=sl->data; + geom_polygon_unique(gs); + struct coord *c=gs->first; while (c <= gs->last) { if (first) { @@ -228,6 +481,14 @@ item_bin_add_coord(ib, gs->first, gs->last-gs->first+1); item_bin_write(ib, f); } + if (ways_split && gs->type==geom_poly_segment_type_way_outer) { + struct item_bin *ib=item_bin_relation_area; + item_bin_add_coord(ib, gs->first, gs->last-gs->first+1); + to_write = 1; + } + if (ways_split && gs->type==geom_poly_segment_type_way_inner) { + have_inner = 1; + } if (boundary->country) { if (!coord_is_equal(*gs->first,*gs->last)) { if (!fu) { @@ -246,16 +507,103 @@ } sl=g_list_next(sl); } + if (ways_split && to_write ) { + struct item_bin *ib=item_bin_relation_area; + if(0 == self_intersect_test(ib)) { //one outer + int i=0; + if(have_inner) { + GList*grouped_segments = + geom_poly_segments_group(boundary->sorted_segments, grouped_segments); + grouped_segments = remove_intersecting_segments(grouped_segments); +#if 1 +printf("OSMID:%d\n",item_bin_get_relationid(boundary->ib)); +fflush(stdout); +#endif + struct attr*triangle_attr = create_triangle_attr(grouped_segments); + + if(triangle_attr && triangle_attr->u.data) { + item_bin_add_attr(ib,triangle_attr); + } + if(triangle_attr) { + g_free(triangle_attr->u.data); + g_free(triangle_attr); + } + } + + i = 0; + while(boundary->area_item_types[i]!=type_none) { //write all matches + ib->type = boundary->area_item_types[i++]; + item_bin_write(ib, ways_split); + }; + } else { /*maybe disjunct outer polygons*/ + GList*grouped_segments = + geom_poly_segments_group(boundary->sorted_segments, grouped_segments); + grouped_segments = remove_intersecting_segments(grouped_segments); + + GList*l = grouped_segments; + int self_intersecting_grouped = 0; + int i=0; + int changed=0; + while (l) { + int changed=0; + struct geom_poly_segment *gs=l->data; + if(gs->type==geom_poly_segment_type_way_outer) { + struct item_bin *ib=item_bin_relation_area; + item_bin_init(ib,boundary->area_item_types[0]); + + if(gs->last-gs->first+1 > 0) { + item_bin_add_coord(ib, gs->first, gs->last-gs->first+1); + } + if( 0!=self_intersect_test(ib)) { + self_intersecting_grouped = 1; + } else { +#if 0 +printf("OSMID2:%d boundary->area_item_types[0]=%X\n",item_bin_get_relationid(boundary->ib),boundary->area_item_types[0]); +fflush(stdout); +#endif + struct attr *triangle_attr = create_triangle_attr(grouped_segments); + + if(triangle_attr && triangle_attr->u.data) { + int i=0; + item_bin_add_attr(ib,triangle_attr); + g_free(triangle_attr->u.data); + } + i=0; + while(boundary->area_item_types[i]!=type_none) { //write all match + ib->type = boundary->area_item_types[i]; + item_bin_write(ib, ways_split); + ++i; + }; + } + //remove first outer + changed=1; + l=g_list_next(l); + grouped_segments = g_list_remove(grouped_segments,gs); + } + if(!changed) { + l=g_list_next(l); + } + } + + if(self_intersecting_grouped) { + osm_warning("relation",item_bin_get_relationid(boundary->ib),0,"Self intersecting relation area\n"); + } + g_list_free(grouped_segments); + } + } + ret=process_boundaries_insert(ret, boundary); l=g_list_next(l); if (f) fclose(f); + if (ways_split) + fclose(ways_split); if (fu) { if (boundary->country) osm_warning("relation",item_bin_get_relationid(boundary->ib),0,"Broken country polygon '%s'\n",boundary->iso2); fclose(fu); } - + g_free(boundary->area_item_types); } #if 0 printf("hierarchy\n"); Index: maptool/CMakeLists.txt =================================================================== --- maptool/CMakeLists.txt (revision 5005) +++ maptool/CMakeLists.txt (working copy) @@ -2,7 +2,7 @@ if(BUILD_MAPTOOL) add_definitions( -DMODULE=maptool ${NAVIT_COMPILE_FLAGS}) include_directories(${CMAKE_CURRENT_SOURCE_DIR}) - 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) + 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) if(NOT MSVC) 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) endif(NOT MSVC) Index: maptool/osm.c =================================================================== --- maptool/osm.c (revision 5005) +++ maptool/osm.c (working copy) @@ -672,6 +672,7 @@ "w natural=scree poly_scree\n" "w natural=scrub poly_scrub\n" "w natural=water poly_water\n" + "w natural=wetland poly_wetland\n" "w natural=wood poly_wood\n" "w piste:type=downhill,piste:difficulty=advanced piste_downhill_advanced\n" "w piste:type=downhill,piste:difficulty=easy piste_downhill_easy\n" @@ -709,6 +710,7 @@ "w waterway=canal water_canal\n" "w waterway=drain water_drain\n" "w waterway=river water_river\n" + "w water=river water_river\n" "w waterway=riverbank poly_water\n" "w waterway=stream water_stream\n" "w barrier=ditch ditch\n" @@ -942,6 +944,24 @@ char buffer[BUFFER_SIZE*2+2]; if (in_relation) { relation_add_tag(k,v); + + //for relation areas we don't set flags + strcpy(buffer,"*=*"); + if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer))) + attr_present[idx]=1; + + sprintf(buffer,"%s=*", k); + if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer))) + attr_present[idx]=2; + + sprintf(buffer,"*=%s", v); + if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer))) + attr_present[idx]=2; + + sprintf(buffer,"%s=%s", k, v); + if ((idx=(int)(long)g_hash_table_lookup(attr_hash, buffer))) + attr_present[idx]=4; + return; } if (! strcmp(k,"ele")) @@ -1482,9 +1502,69 @@ } +static int +attr_longest_match(struct attr_mapping **mapping, int mapping_count, enum item_type *types, int types_count) +{ + int i,j,longest=0,ret=0,sum,val; + struct attr_mapping *curr; + for (i = 0 ; i < mapping_count ; i++) { + sum=0; + curr=mapping[i]; + for (j = 0 ; j < curr->attr_present_idx_count ; j++) { + val=attr_present[curr->attr_present_idx[j]]; + if (val) + sum+=val; + else { + sum=-1; + break; + } + } + if (sum > longest) { + longest=sum; + ret=0; + } + if (sum > 0 && sum == longest && ret < types_count) + types[ret++]=curr->type; + } + return ret; +} + +static void +attr_longest_match_clear(void) +{ + memset(attr_present, 0, sizeof(*attr_present)*attr_present_count); +} + void osm_end_relation(struct maptool_osm *osm) { + int count,itcount=0; + enum item_type types[10]; + types[0] = type_none; + + count=attr_longest_match(attr_mapping_way, attr_mapping_way_count, types, sizeof(types)/sizeof(enum item_type)); + if(count>0) { + int ii; + + struct attr attr; + attr.type = attr_item_types; + attr.u.item_types = NULL; + + for(ii=0;ii0) { + item_bin_add_attr(item_bin, &attr); + } + } + attr_longest_match_clear(); + in_relation=0; if ((!strcmp(relation_type, "multipolygon") || !strcmp(relation_type, "boundary")) && boundary) { #if 0 @@ -1553,40 +1633,6 @@ } } - -static int -attr_longest_match(struct attr_mapping **mapping, int mapping_count, enum item_type *types, int types_count) -{ - int i,j,longest=0,ret=0,sum,val; - struct attr_mapping *curr; - for (i = 0 ; i < mapping_count ; i++) { - sum=0; - curr=mapping[i]; - for (j = 0 ; j < curr->attr_present_idx_count ; j++) { - val=attr_present[curr->attr_present_idx[j]]; - if (val) - sum+=val; - else { - sum=-1; - break; - } - } - if (sum > longest) { - longest=sum; - ret=0; - } - if (sum > 0 && sum == longest && ret < types_count) - types[ret++]=curr->type; - } - return ret; -} - -static void -attr_longest_match_clear(void) -{ - memset(attr_present, 0, sizeof(*attr_present)*attr_present_count); -} - void osm_end_way(struct maptool_osm *osm) { @@ -2562,6 +2608,7 @@ processed_nodes=processed_nodes_out=processed_ways=processed_relations=processed_tiles=0; sig_alrm(0); while ((ib=read_item(in))) { + int self_intersect; #if 0 fprintf(stderr,"type 0x%x len %d clen %d\n", ib->type, ib->len, ib->clen); #endif @@ -2577,7 +2624,8 @@ if (ni) { c[i]=ni->c; if (ni->ref_way > 1 && i != 0 && i != ccount-1 && i != last && item_get_default_flags(ib->type)) { - write_item_part(out, out_index, out_graph, ib, last, i, &last_id); + if(!(0 != self_intersect_test(ib) && item_type_is_area(ib->type))) + write_item_part(out, out_index, out_graph, ib, last, i, &last_id); last=i; } } else if (final) { @@ -2592,10 +2640,16 @@ } } } + self_intersect = self_intersect_test(ib); + if(self_intersect!=0 && (/*ib->type == type_water_line ||*/ item_type_is_area(ib->type)) ) { + osm_warning("way", item_bin_get_wayid(ib), 0, "Self intersecting area\n"); + } if (ccount) { - write_item_part(out, out_index, out_graph, ib, last, ccount-1, &last_id); + if(!(0 != self_intersect && item_type_is_area(ib->type))) + write_item_part(out, out_index, out_graph, ib, last, ccount-1, &last_id); if (final && ib->type == type_water_line && out_coastline) { - write_item_part(out_coastline, NULL, NULL, ib, last, ccount-1, NULL); + //if(0 == self_intersect) + write_item_part(out_coastline, NULL, NULL, ib, last, ccount-1, NULL); } } } Index: maptool/maptool.c =================================================================== --- maptool/maptool.c (revision 5005) +++ maptool/maptool.c (working copy) @@ -745,7 +745,6 @@ } } - int main(int argc, char **argv) { #if 0 @@ -821,7 +820,7 @@ return 0; } phase=0; - +/**/ // input from an OSM file if (p.input == 0) { if (start_phase(&p, "collecting data")) { @@ -857,6 +856,7 @@ if (start_phase(&p,"generating coastlines")) { osm_process_coastlines(&p, suffix); } +/**/ if (start_phase(&p,"assinging towns to countries")) { FILE *towns=tempfile(suffix,"towns",0),*boundaries=NULL,*ways=NULL; if (towns) { Index: attr_def.h =================================================================== --- attr_def.h (revision 5005) +++ attr_def.h (working copy) @@ -430,6 +430,7 @@ ATTR(announcement) ATTR(cursor) ATTR(config) +ATTR(triangles) ATTR2(0x0008ffff,type_object_end) ATTR2(0x00090000,type_coord_begin) ATTR2(0x0009ffff,type_coord_end) Index: attr.c =================================================================== --- attr.c (revision 5005) +++ attr.c (working copy) @@ -492,6 +492,10 @@ int attr_data_size(struct attr *attr) { + if (attr->type == attr_triangles) { + int*pi=(u32*)(attr->u.data); + return sizeof(u32)*(2 + 3*pi[1] + pi[0]); + } if (attr->type >= attr_type_string_begin && attr->type <= attr_type_string_end) return strlen(attr->u.str)+1; if (attr->type >= attr_type_int_begin && attr->type <= attr_type_int_end) Index: graphics.c =================================================================== --- graphics.c (revision 5005) +++ graphics.c (working copy) @@ -50,6 +50,7 @@ #include "callback.h" #include "file.h" #include "event.h" +#include "types.h" //############################################################################################################## @@ -452,8 +453,6 @@ */ void graphics_gc_destroy(struct graphics_gc *gc) { - if (!gc) - return; gc->meth.gc_destroy(gc->priv); g_free(gc); } @@ -860,6 +859,7 @@ char *label; int z_order; int count; + struct attr *triangles;/*stores triangle data for polygons that need triangulation (polygons with holes)*/ struct coord c[0]; }; @@ -876,6 +876,7 @@ struct displayitem *di=dl->hash_entries[i].di; while (di) { struct displayitem *next=di->next; + g_free(di->triangles); g_free(di); di=next; } @@ -889,7 +890,7 @@ * @returns <> * @author Martin Schaller (04/2008) */ -static void display_add(struct hash_entry *entry, struct item *item, int count, struct coord *c, char **label, int label_count) +static void display_add(struct hash_entry *entry, struct item *item, int count, struct coord *c, char **label, int label_count, struct attr*poly_attr) { struct displayitem *di; int len,i; @@ -910,6 +911,7 @@ p+=sizeof(*di)+count*sizeof(*c); di->item=*item; di->z_order=0; + di->triangles = poly_attr; if (label && label_count) { di->label=p; for (i = 0 ; i < label_count ; i++) { @@ -1727,6 +1729,7 @@ char *path; while (di) { + int i,count=di->count,mindist=dc->mindist; di->z_order=++(gra->current_z_order); @@ -1741,12 +1744,50 @@ if (dc->type == type_poly_water_tiled) mindist=0; if (dc->e->type == element_polyline) - count=transform(dc->trans, dc->pro, di->c, pa, count, mindist, e->u.polyline.width, width); + count=transform(dc->trans, dc->pro, di->c, pa, count, 0/*mindist*/, e->u.polyline.width, width); else - count=transform(dc->trans, dc->pro, di->c, pa, count, mindist, 0, NULL); + count=transform(dc->trans, dc->pro, di->c, pa, count, 0/*mindist*/, 0, NULL); switch (e->type) { case element_polygon: - graphics_draw_polygon_clipped(gra, gc, pa, count); + if(di->triangles) { + u32*pi=di->triangles->u.data; + int ii; +#if 0 +int jj; +for (jj=0;jjtrans, dc->pro, ca, pa, 3, 0, 0, NULL); + //TODO use some caching mechanism to avoid unnecessary recaculations + /* + if( + (pa[0].xr.lu.x || gra->r.rl.xr.lu.y || gra->r.rl.yr.lu.x || gra->r.rl.xr.lu.y || gra->r.rl.yr.lu.x || gra->r.rl.xr.lu.y || gra->r.rl.ymr) { while ((item=map_rect_get_item(displaylist->mr))) { + + + + + + + int label_count=0; char *labels[2]; struct hash_entry *entry; @@ -2110,6 +2158,22 @@ dbg(0,"point count overflow %d for %s "ITEM_ID_FMT"\n", count,item_to_name(item->type),ITEM_ID_ARGS(*item)); displaylist->dc.maxlen=max*2; } + + struct attr* poly_attr = g_new0(struct attr,1); + + if(item_attr_get(item,attr_triangles,poly_attr)) { + int *pi = poly_attr->u.data; + int *pi2; + int ii; + pi2 = g_malloc(attr_data_size(poly_attr)); + memcpy(pi2,pi,attr_data_size(poly_attr)); + poly_attr->u.data = pi2; + + } else { + g_free(poly_attr); + poly_attr = NULL; + } + if (item_is_custom_poi(*item)) { if (item_attr_get(item, attr_icon_src, &attr2)) labels[1]=map_convert_string(displaylist->m, attr2.u.str); @@ -2128,10 +2192,10 @@ labels[0]=NULL; if (displaylist->conv && label_count) { labels[0]=map_convert_string(displaylist->m, labels[0]); - display_add(entry, item, count, ca, labels, label_count); + display_add(entry, item, count, ca, labels, label_count ,poly_attr); map_convert_free(labels[0]); } else - display_add(entry, item, count, ca, labels, label_count); + display_add(entry, item, count, ca, labels, label_count ,poly_attr); if (labels[1]) map_convert_free(labels[1]); workload++;