From 27c31404caa399f424d742e1619b4a5426999cad Mon Sep 17 00:00:00 2001 From: PSadrozinski <psadrozinski@gmail.com> Date: Sun, 25 Dec 2011 20:54:37 -0500 Subject: [PATCH] Compile warnings cleanup round 1 --- src/Airports/GenAirports/build.cxx | 12 +- src/Lib/Geometry/util.cxx | 18 +- src/Lib/Polygon/clipper.cpp | 2 + src/Lib/e00/e00.cxx | 19 +- src/Lib/poly2tri/construct.c | 1597 +++++++++++++--------------- src/Lib/poly2tri/monotone.c | 2 +- src/Lib/shapelib/dbfadd.c | 4 +- src/Lib/shapelib/shputils.c | 335 +++--- src/Lib/vpf/table.cxx | 64 +- src/Lib/vpf/value.cxx | 62 +- src/Prep/GDALChop/gdalchop.cxx | 630 +++++------ src/Prep/Terra/Subdivision.cc | 31 +- 12 files changed, 1369 insertions(+), 1407 deletions(-) diff --git a/src/Airports/GenAirports/build.cxx b/src/Airports/GenAirports/build.cxx index 1d257ce5..9e074dc5 100644 --- a/src/Airports/GenAirports/build.cxx +++ b/src/Airports/GenAirports/build.cxx @@ -605,7 +605,7 @@ void build_airport( string airport_id, float alt_m, } #else /* Ralf Gerlich: Generate Taxiways in specified order from bottom to top */ - for ( i=0; i<taxiways.size(); ++i ) { + for ( i=0; i<(int)taxiways.size(); ++i ) { SG_LOG( SG_GENERAL, SG_DEBUG, "generating " << i ); build_runway( taxiways[i], alt_m, &rwy_polys, &texparams, &accum, @@ -970,7 +970,7 @@ void build_airport( string airport_id, float alt_m, tri_materials.push_back( "Grass" ); std::vector < SGGeod > geodNodes; - for ( j = 0; j < nodes.get_node_list().size(); j++ ) { + for ( j = 0; j < (int)nodes.get_node_list().size(); j++ ) { Point3D node = nodes.get_node_list()[j]; geodNodes.push_back( SGGeod::fromDegM( node.x(), node.y(), node.z() ) ); } @@ -1162,7 +1162,7 @@ void build_airport( string airport_id, float alt_m, strip_materials.push_back( "Grass" ); std::vector < SGGeod > geodNodes; - for ( j = 0; j < nodes.get_node_list().size(); j++ ) { + for ( j = 0; j < (int)nodes.get_node_list().size(); j++ ) { Point3D node = nodes.get_node_list()[j]; geodNodes.push_back( SGGeod::fromDegM( node.x(), node.y(), node.z() ) ); } @@ -1285,7 +1285,7 @@ void build_airport( string airport_id, float alt_m, wgs84_nodes.push_back( cart ); } SGSphered d; - for ( i = 0; i < wgs84_nodes.size(); ++i ) { + for ( i = 0; i < (int)wgs84_nodes.size(); ++i ) { d.expandBy(wgs84_nodes[ i ]); } @@ -1306,13 +1306,13 @@ void build_airport( string airport_id, float alt_m, string name = airport_id + ".btg"; std::vector< SGVec3f > normals_3f; - for ( i=0; i < normals.get_node_list().size(); i++ ) { + for ( i=0; i < (int)normals.get_node_list().size(); i++ ) { Point3D node = normals.get_node_list()[i]; normals_3f.push_back( node.toSGVec3f() ); } std::vector< SGVec2f > texcoords_2f; - for ( i=0; i < texcoords.get_node_list().size(); i++ ) { + for ( i=0; i < (int)texcoords.get_node_list().size(); i++ ) { Point3D node = texcoords.get_node_list()[i]; texcoords_2f.push_back( node.toSGVec2f() ); } diff --git a/src/Lib/Geometry/util.cxx b/src/Lib/Geometry/util.cxx index 6e335da1..0c7f3427 100644 --- a/src/Lib/Geometry/util.cxx +++ b/src/Lib/Geometry/util.cxx @@ -69,7 +69,7 @@ getIntersection (const Point3D &p0, const Point3D &p1, void makePolygon (const Point3D &p, double width, TGPolygon &polygon) { - double x, y, az; + double x = 0.0f, y = 0.0f, az = 0.0f; double lon = p.x(); double lat = p.y(); @@ -101,7 +101,12 @@ makePolygon (const Line &line, double width, TGPolygon &polygon) const Point3D p1 = line.getPoint(i); const Point3D p2 = line.getPoint(i+1); - double angle1, angle2, dist, x, y, az; + double angle1 = 0.0f; + double angle2 = 0.0f; + double dist = 0.0f; + double x = 0.0f; + double y = 0.0f; + double az = 0.0f; geo_inverse_wgs_84(0, p1.y(), p1.x(), p2.y(), p2.x(), &angle1, &angle2, &dist); polygon.erase(); @@ -420,9 +425,12 @@ makePolygonsTP (const Line &line, double width, poly_list& polys, texparams_list Point3D prev_inner = Point3D(0.0f, 0.0f, 0.0f); Point3D prev_outer = Point3D(0.0f, 0.0f, 0.0f); - double last_end_v; - double heading, az2, dist; - double pt_x, pt_y; + double last_end_v = 0.0f; + double heading = 0.0f; + double az2 = 0.0f; + double dist = 0.0f; + double pt_x = 0.0f; + double pt_y = 0.0f; TGPolygon poly; TGTexParams tp; diff --git a/src/Lib/Polygon/clipper.cpp b/src/Lib/Polygon/clipper.cpp index 3cc2e58a..da08072a 100644 --- a/src/Lib/Polygon/clipper.cpp +++ b/src/Lib/Polygon/clipper.cpp @@ -415,6 +415,8 @@ double Area(const Polygon &poly) case rtHi: UseFullInt64Range = true; break; + case rtLo: + break; case rtError: throw "Coordinate exceeds range bounds."; } diff --git a/src/Lib/e00/e00.cxx b/src/Lib/e00/e00.cxx index b2a7247d..4ebd385d 100644 --- a/src/Lib/e00/e00.cxx +++ b/src/Lib/e00/e00.cxx @@ -37,7 +37,7 @@ strAppend (string &s, int i) s += buf; } - +#if 0 /* UNUSED */ /** * Append a double-precision real to a string. */ @@ -48,7 +48,7 @@ strAppend (string &s, double f) sprintf(buf, "%f", f); s += buf; } - +#endif /** * Skip newlines. @@ -141,7 +141,6 @@ checkZeros (istream &input) } } - /** * Check that the expected integer appears, or throw an exception. */ @@ -157,7 +156,7 @@ expect (istream &input, int i) } } - +#if 0 /* UNUSED */ /** * Check that the expected real number appears, or throw an exception. */ @@ -172,7 +171,7 @@ expect (istream &input, double f) throw E00Exception(message.c_str()); } } - +#endif /** * Check that the expected string appears, or throw an exception. @@ -541,8 +540,6 @@ E00::readIFO () { int line_pos = 0; string line = ""; - int intval; - double realval; checkPrecision(*_input); @@ -689,7 +686,7 @@ E00::postProcess () const E00::IFO * E00::getIFO (const string &fileName) const { - for (int i = 0; i < ifo_section.size(); i++) { + for (int i = 0; i < (int)ifo_section.size(); i++) { if (ifo_section[i].fileName == fileName) return &(ifo_section[i]); } @@ -705,7 +702,7 @@ E00::getIFOItem (const string &fileName, int entry, return 0; int pos = -1; - for (int i = 0; i < ifo->defs.size(); i++) { + for (int i = 0; i < (int)ifo->defs.size(); i++) { if (ifo->defs[i].itemName == itemName) pos = i; } @@ -723,8 +720,8 @@ E00::getIFOItemType (const string &fileName, const string &itemName) const if (ifo == 0) return 0; - int pos = -1; - for (int i = 0; i < ifo->defs.size(); i++) { + // int pos = -1; + for (int i = 0; i < (int)ifo->defs.size(); i++) { if (ifo->defs[i].itemName == itemName) return &(ifo->defs[i].itemType); } diff --git a/src/Lib/poly2tri/construct.c b/src/Lib/poly2tri/construct.c index 24745a19..20a499ff 100644 --- a/src/Lib/poly2tri/construct.c +++ b/src/Lib/poly2tri/construct.c @@ -6,9 +6,9 @@ # include <memory.h> #endif -node_t qs[QSIZE]; /* Query structure */ -trap_t tr[TRSIZE]; /* Trapezoid structure */ -segment_t seg[SEGSIZE]; /* Segment table */ +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; @@ -16,209 +16,201 @@ 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; + 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; + 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; +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; + 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; +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; + 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; +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); + 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; +point_t * v0; +point_t *v1; { - return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x)); + 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; +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); + 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; +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); + 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) +/* 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 - * \ + * \ + * 1 \ 2 + * \ * ----------------------------------- * 3 */ static int init_query_structure(segnum) - int segnum; +int segnum; { - int i1, i2, i3, i4, i5, i6, i7, root; - int t1, t2, t3, t4; - segment_t *s = &seg[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)); + 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; + 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].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[i1].left = i3 = newnode(); + qs[i3].nodetype = T_Y; + _min(&qs[i3].yval, &s->v0, &s->v1); /* root */ + qs[i3].parent = i1; - qs[i5].right = i7 = newnode(); - qs[i7].nodetype = T_SINK; - qs[i7].parent = i5; + qs[i3].left = i4 = newnode(); + qs[i4].nodetype = T_SINK; + qs[i4].parent = i3; - t1 = newtrap(); /* middle left */ - t2 = newtrap(); /* middle right */ - t3 = newtrap(); /* bottom-most */ - t4 = newtrap(); /* topmost */ + qs[i3].right = i5 = newnode(); + qs[i5].nodetype = T_X; + qs[i5].segnum = segnum; + qs[i5].parent = i3; - 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; + qs[i5].left = i6 = newnode(); + qs[i6].nodetype = T_SINK; + qs[i6].parent = i5; - tr[t1].state = tr[t2].state = ST_VALID; - tr[t3].state = tr[t4].state = ST_VALID; + qs[i5].right = i7 = newnode(); + qs[i7].nodetype = T_SINK; + qs[i7].parent = i5; - qs[i2].trnum = t4; - qs[i4].trnum = t3; - qs[i6].trnum = t1; - qs[i7].trnum = t2; + t1 = newtrap(); /* middle left */ + t2 = newtrap(); /* middle right */ + t3 = newtrap(); /* bottom-most */ + t4 = newtrap(); /* topmost */ - s->is_inserted = TRUE; - return root; + 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; } @@ -228,55 +220,44 @@ static int init_query_structure(segnum) */ static int is_left_of(segnum, v) - int segnum; - point_t *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)); + 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)); } - 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; + + if (area > 0.0) + return TRUE; + else + return FALSE; } @@ -286,74 +267,69 @@ static int is_left_of(segnum, v) /* whether the segment which shares this endpoint is already inserted */ static int inserted(segnum, whichpt) - int segnum; - int whichpt; +int segnum; +int whichpt; { - if (whichpt == FIRSTPT) - return seg[seg[segnum].prev].is_inserted; - else - return seg[seg[segnum].next].is_inserted; + 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. +/* 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; +point_t * v; +point_t *vo; +int r; { - node_t *rptr = &qs[r]; - - switch (rptr->nodetype) - { + node_t *rptr = &qs[r]; + + switch (rptr->nodetype) { case T_SINK: - return rptr->trnum; - + 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 */ + 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 */ + 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; + fprintf(stderr, "Haggu !!!!!\n"); + break; } + + return -1; } -/* Thread in the segment into the existing trapezoidation. The +/* 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 @@ -361,67 +337,65 @@ int locate_endpoint(v, vo, r) */ static int merge_trapezoids(segnum, tfirst, tlast, side) - int segnum; - int tfirst; - int tlast; - int side; +int segnum; +int tfirst; +int tlast; +int side; { - int t, tnext, cond; - int ptnext; + 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; - /* 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; + + return 0; } @@ -432,575 +406,483 @@ static int merge_trapezoids(segnum, tfirst, tlast, side) */ static int add_segment(segnum) - int segnum; +int segnum; { - segment_t s; - segment_t *so = &seg[segnum]; - int tu, tl, sk, tfirst, tlast; - int tfirstr, tlastr, tfirstl, tlastl; - int i1, i2, t, tn; - point_t tpt; - int tritop = 0, tribot = 0, is_swapped = 0; - int tmptriseg; + segment_t s; + int tu, tl, sk, tfirst, tlast; + int tfirstr = -1, tlastr = -1, tfirstl = -1, tlastl = -1; + int i1, i2, t, 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; + 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; + 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; + 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].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; + 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; + /* Now update the query structure and obtain the sinks for the */ + /* two trapezoids */ - qs[i1].nodetype = T_SINK; - qs[i1].trnum = tu; - qs[i1].parent = sk; + i1 = newnode(); /* Upper trapezoid sink */ + i2 = newnode(); /* Lower trapezoid sink */ + sk = tr[tu].sink; - qs[i2].nodetype = T_SINK; - qs[i2].trnum = tl; - qs[i2].parent = sk; + 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; - 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; + 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; + 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); + 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; + 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].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; + 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; - qs[i1].nodetype = T_SINK; - qs[i1].trnum = tu; - qs[i1].parent = sk; + /* Now update the query structure and obtain the sinks for the */ + /* two trapezoids */ - qs[i2].nodetype = T_SINK; - qs[i2].trnum = tl; - qs[i2].parent = sk; + i1 = newnode(); /* Upper trapezoid sink */ + i2 = newnode(); /* Lower trapezoid sink */ + sk = tr[tu].sink; - tr[tu].sink = i1; - tr[tl].sink = i2; - tlast = tu; + 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; } - 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; + /* Thread the segment into the query tree creating a new X-node */ + /* First, split all the trapezoids which are intersected by s into */ + /* two */ - qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */ - qs[i1].trnum = t; - qs[i1].parent = sk; + t = tfirst; /* topmost trapezoid */ - qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */ - qs[i2].trnum = tn = newtrap(); - tr[tn].state = ST_VALID; - qs[i2].parent = sk; + 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 */ - if (t == tfirst) - tfirstr = tn; - if (_equal_to(&tr[t].lo, &tr[tlast].lo)) - tlastr = tn; + qs[sk].nodetype = T_X; + qs[sk].segnum = segnum; + qs[sk].left = i1; + qs[sk].right = i2; - tr[tn] = tr[t]; - tr[t].sink = i1; - tr[tn].sink = i2; - t_sav = t; - tn_sav = tn; + qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */ + qs[i1].trnum = t; + qs[i1].parent = sk; - /* error */ + qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */ + qs[i2].trnum = tn = newtrap(); + tr[tn].state = ST_VALID; + qs[i2].parent = sk; - if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */ - { - fprintf(stderr, "add_segment: error\n"); - 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; + if (t == tfirst) + tfirstr = tn; + if (_equal_to(&tr[t].lo, &tr[tlast].lo)) + tlastr = tn; - 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 */ + tr[tn] = tr[t]; + tr[t].sink = i1; + tr[tn].sink = i2; + t_sav = t; + tn_sav = tn; - 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; - } + /* error */ + if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) { /* case cannot arise */ + fprintf(stderr, "add_segment: error\n"); + 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; - 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[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; + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } - if (is_swapped) - tmptriseg = seg[segnum].prev; - else - tmptriseg = seg[segnum].next; + 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 ((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; - } + 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 */ - /* 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; + if (is_swapped) + tmptriseg = seg[segnum].prev; + else + tmptriseg = seg[segnum].next; - 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 ((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; + } - 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; + 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[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[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[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; + tr[tr[t].u0].d0 = t; + tr[tr[t].u1].d0 = t; + tr[tr[tn].u0].d0 = tn; + } - 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; + 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; + } + } - /* 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; + 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 tmptriseg; + + 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].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{ + 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); + /* 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 */ - seg[segnum].is_inserted = TRUE; - return 0; + 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; } @@ -1009,51 +891,50 @@ static int add_segment(segnum) * the segment is inserted into the trapezoidation subsequently */ static int find_new_roots(segnum) - int segnum; +int segnum; { - segment_t *s = &seg[segnum]; - - if (s->is_inserted) + 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; - - 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; +int nseg; { - register int i; - int root, h; - - /* Add the first segment and get the query structure and trapezoid */ - /* list initialised */ + register int i; + int root, h; - root = init_query_structure(choose_segment()); + /* Add the first segment and get the query structure and trapezoid */ + /* list initialised */ - 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++) - add_segment(choose_segment()); - - /* Find a new root for each of the segment endpoints */ - for (i = 1; i <= nseg; i++) - find_new_roots(i); + 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++) + add_segment(choose_segment()); + + /* 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++) - add_segment(choose_segment()); - return 0; + for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++) + add_segment(choose_segment()); + + return 0; } diff --git a/src/Lib/poly2tri/monotone.c b/src/Lib/poly2tri/monotone.c index 9552302f..7138ee3e 100644 --- a/src/Lib/poly2tri/monotone.c +++ b/src/Lib/poly2tri/monotone.c @@ -289,7 +289,7 @@ static int traverse_polygon(mcur, trnum, from, dir) trap_t *t = &tr[trnum]; int mnew; int v0, v1; - int retval; + int retval = 0; int do_switch = FALSE; if ((trnum <= 0) || visited[trnum]) diff --git a/src/Lib/shapelib/dbfadd.c b/src/Lib/shapelib/dbfadd.c index f9782087..c62e8ee2 100644 --- a/src/Lib/shapelib/dbfadd.c +++ b/src/Lib/shapelib/dbfadd.c @@ -57,9 +57,11 @@ static char rcsid[] = "$Id: dbfadd.c,v 1.7 2002/01/15 14:36:07 warmerda Exp $"; -#include "shapefil.h" #include <math.h> #include <stdlib.h> +#include <stdio.h> + +#include "shapefil.h" int main( int argc, char ** argv ) diff --git a/src/Lib/shapelib/shputils.c b/src/Lib/shapelib/shputils.c index 38f59c88..7372feef 100644 --- a/src/Lib/shapelib/shputils.c +++ b/src/Lib/shapelib/shputils.c @@ -76,8 +76,10 @@ * Initial revision */ +#if 0 static char rcsid[] = "$Id: shputils.c,v 1.7 2003/02/25 17:20:22 warmerda Exp $"; +#endif #include <stdlib.h> #include "shapefil.h" @@ -121,6 +123,8 @@ void setext(char *pt, char *ext); int strncasecmp2(char *s1, char *s2, int n); void mergefields(void); void findselect(void); +int findunit(char *unit); + void showitems(void); int selectrec(); int check_theme_bnd(); @@ -394,6 +398,9 @@ int main( int argc, char ** argv ) DBFWriteDoubleAttribute(hDBFappend, jRecord, pt[i], (DBFReadDoubleAttribute( hDBF, iRecord, i )) ); break; + + default: + break; } } } @@ -686,8 +693,8 @@ char *pt; isum = isum + itmp; } mean=isum/maxrec; - if (ilow < ihigh) printf("%d to %d \t(%.1f)",ilow,ihigh,mean); - else if (ilow == ihigh) printf("= %d",ilow); + if (ilow < ihigh) printf("%ld to %ld \t(%.1f)",ilow,ihigh,mean); + else if (ilow == ihigh) printf("= %ld",ilow); else printf("No Values"); break; @@ -712,6 +719,8 @@ char *pt; else printf("No Values"); break; + default: + break; } } @@ -720,161 +729,207 @@ char *pt; int selectrec() { -long int value, ty; + long int value, ty; - ty = DBFGetFieldInfo( hDBF, iselectitem, NULL, &iWidth, &iDecimals); - switch(ty) - { - case FTString: - puts("Invalid Item"); - iselect=FALSE; - break; - case FTInteger: - value = DBFReadIntegerAttribute( hDBF, iRecord, iselectitem ); - for (j = 0; j<selcount; j++) - { - if (selectvalues[j] == value) - if (iunselect) return(0); /* Keep this record */ - else return(1); /* Skip this record */ - } - break; - case FTDouble: - puts("Invalid Item"); - iselect=FALSE; - break; - } - if (iunselect) return(1); /* Skip this record */ - else return(0); /* Keep this record */ + ty = DBFGetFieldInfo( hDBF, iselectitem, NULL, &iWidth, &iDecimals); + switch(ty) { + case FTString: + puts("Invalid Item"); + iselect=FALSE; + break; + + case FTInteger: + value = DBFReadIntegerAttribute( hDBF, iRecord, iselectitem ); + for (j = 0; j<selcount; j++) { + if (selectvalues[j] == value) { + if (iunselect) { + return(0); /* Keep this record */ + } else { + return(1); /* Skip this record */ + } + } + } + break; + + case FTDouble: + puts("Invalid Item"); + iselect=FALSE; + break; + + default: + break; + } + + if (iunselect) { + return(1); /* Skip this record */ + } else { + return(0); /* Keep this record */ + } } int check_theme_bnd() { if ( (adfBoundsMin[0] >= cxmin) && (adfBoundsMax[0] <= cxmax) && - (adfBoundsMin[1] >= cymin) && (adfBoundsMax[1] <= cymax) ) - { /** Theme is totally inside clip area **/ - if (ierase) nEntities=0; /** SKIP THEME **/ - else iclip=FALSE; /** WRITE THEME (Clip not needed) **/ + (adfBoundsMin[1] >= cymin) && (adfBoundsMax[1] <= cymax) ) { + /** Theme is totally inside clip area **/ + if (ierase) { + nEntities=0; /** SKIP THEME **/ + } else { + iclip=FALSE; /** WRITE THEME (Clip not needed) **/ + } } if ( ( (adfBoundsMin[0] < cxmin) && (adfBoundsMax[0] < cxmin) ) || ( (adfBoundsMin[1] < cymin) && (adfBoundsMax[1] < cymin) ) || ( (adfBoundsMin[0] > cxmax) && (adfBoundsMax[0] > cxmax) ) || - ( (adfBoundsMin[1] > cymax) && (adfBoundsMax[1] > cymax) ) ) - { /** Theme is totally outside clip area **/ - if (ierase) iclip=FALSE; /** WRITE THEME (Clip not needed) **/ - else nEntities=0; /** SKIP THEME **/ + ( (adfBoundsMin[1] > cymax) && (adfBoundsMax[1] > cymax) ) ) { + /** Theme is totally outside clip area **/ + if (ierase) { + iclip=FALSE; /** WRITE THEME (Clip not needed) **/ + } else { + nEntities=0; /** SKIP THEME **/ + } } - if (nEntities == 0) + if (nEntities == 0) { puts("WARNING: Theme is outside the clip area."); /** SKIP THEME **/ + } + + return 0; } -clip_boundary() +int clip_boundary() { int inside; int prev_outside; int i2; int j2; - /*** FIRST check the boundary of the feature ***/ - if ( ( (psCShape->dfXMin < cxmin) && (psCShape->dfXMax < cxmin) ) || - ( (psCShape->dfYMin < cymin) && (psCShape->dfYMax < cymin) ) || - ( (psCShape->dfXMin > cxmax) && (psCShape->dfXMax > cxmax) ) || - ( (psCShape->dfYMin > cymax) && (psCShape->dfYMax > cymax) ) ) - { /** Feature is totally outside clip area **/ - if (ierase) return(1); /** WRITE RECORD **/ - else return(0); /** SKIP RECORD **/ - } + /*** FIRST check the boundary of the feature ***/ + if ( ( (psCShape->dfXMin < cxmin) && (psCShape->dfXMax < cxmin) ) || + ( (psCShape->dfYMin < cymin) && (psCShape->dfYMax < cymin) ) || + ( (psCShape->dfXMin > cxmax) && (psCShape->dfXMax > cxmax) ) || + ( (psCShape->dfYMin > cymax) && (psCShape->dfYMax > cymax) ) ) { + /** Feature is totally outside clip area **/ + if (ierase) { + return(1); /** WRITE RECORD **/ + } else { + return(0); /** SKIP RECORD **/ + } + } - if ( (psCShape->dfXMin >= cxmin) && (psCShape->dfXMax <= cxmax) && - (psCShape->dfYMin >= cymin) && (psCShape->dfYMax <= cymax) ) - { /** Feature is totally inside clip area **/ - if (ierase) return(0); /** SKIP RECORD **/ - else return(1); /** WRITE RECORD **/ - } + if ( (psCShape->dfXMin >= cxmin) && (psCShape->dfXMax <= cxmax) && + (psCShape->dfYMin >= cymin) && (psCShape->dfYMax <= cymax) ) { + /** Feature is totally inside clip area **/ + if (ierase) { + return(0); /** SKIP RECORD **/ + } else { + return(1); /** WRITE RECORD **/ + } + } - if (iinside) - { /** INSIDE * Feature might touch the boundary or could be outside **/ - if (ierase) return(1); /** WRITE RECORD **/ - else return(0); /** SKIP RECORD **/ - } - - if (itouch) - { /** TOUCH **/ - if ( ( (psCShape->dfXMin <= cxmin) || (psCShape->dfXMax >= cxmax) ) && - (psCShape->dfYMin >= cymin) && (psCShape->dfYMax <= cymax) ) - { /** Feature intersects the clip boundary only on the X axis **/ - if (ierase) return(0); /** SKIP RECORD **/ - else return(1); /** WRITE RECORD **/ - } + if (iinside) { /** INSIDE * Feature might touch the boundary or could be outside **/ + if (ierase) { + return(1); /** WRITE RECORD **/ + } else { + return(0); /** SKIP RECORD **/ + } + } - if ( (psCShape->dfXMin >= cxmin) && (psCShape->dfXMax <= cxmax) && - ( (psCShape->dfYMin <= cymin) || (psCShape->dfYMax >= cymax) ) ) - { /** Feature intersects the clip boundary only on the Y axis **/ - if (ierase) return(0); /** SKIP RECORD **/ - else return(1); /** WRITE RECORD **/ - } + if (itouch) { + /** TOUCH **/ + if ( ( (psCShape->dfXMin <= cxmin) || (psCShape->dfXMax >= cxmax) ) && + (psCShape->dfYMin >= cymin) && (psCShape->dfYMax <= cymax) ) { + /** Feature intersects the clip boundary only on the X axis **/ + if (ierase) { + return(0); /** SKIP RECORD **/ + } else { + return(1); /** WRITE RECORD **/ + } + } + + if ( (psCShape->dfXMin >= cxmin) && (psCShape->dfXMax <= cxmax) && + ( (psCShape->dfYMin <= cymin) || (psCShape->dfYMax >= cymax) ) ) { + /** Feature intersects the clip boundary only on the Y axis **/ + if (ierase) { + return(0); /** SKIP RECORD **/ + } else { + return(1); /** WRITE RECORD **/ + } + } - for( j2 = 0; j2 < psCShape->nVertices; j2++ ) - { /** At least one vertex must be inside the clip boundary **/ - if ( (psCShape->padfX[j2] >= cxmin && psCShape->padfX[j2] <= cxmax) || - (psCShape->padfY[j2] >= cymin && psCShape->padfY[j2] <= cymax) ) - if (ierase) return(0); /** SKIP RECORD **/ - else return(1); /** WRITE RECORD **/ - } + for( j2 = 0; j2 < psCShape->nVertices; j2++ ) { + /** At least one vertex must be inside the clip boundary **/ + if ( (psCShape->padfX[j2] >= cxmin && psCShape->padfX[j2] <= cxmax) || + (psCShape->padfY[j2] >= cymin && psCShape->padfY[j2] <= cymax) ) { + if (ierase) { + return(0); /** SKIP RECORD **/ + } else { + return(1); /** WRITE RECORD **/ + } + } + } - /** All vertices are outside the clip boundary **/ - if (ierase) return(1); /** WRITE RECORD **/ - else return(0); /** SKIP RECORD **/ - } /** End TOUCH **/ + /** All vertices are outside the clip boundary **/ + if (ierase) { + return(1); /** WRITE RECORD **/ + } else { + return(0); /** SKIP RECORD **/ + } + } /** End TOUCH **/ - if (icut) - { /** CUT **/ - /*** Check each vertex in the feature with the Boundary and "CUT" ***/ - /*** THIS CODE WAS NOT COMPLETED! READ NOTE AT THE BOTTOM ***/ - i2=0; - prev_outside=FALSE; - for( j2 = 0; j2 < psCShape->nVertices; j2++ ) - { - inside = psCShape->padfX[j2] >= cxmin && psCShape->padfX[j2] <= cxmax && - psCShape->padfY[j2] >= cymin && psCShape->padfY[j2] <= cymax ; + if (icut) { + /** CUT **/ + /*** Check each vertex in the feature with the Boundary and "CUT" ***/ + /*** THIS CODE WAS NOT COMPLETED! READ NOTE AT THE BOTTOM ***/ + i2=0; + prev_outside=FALSE; + for( j2 = 0; j2 < psCShape->nVertices; j2++ ) { + inside = psCShape->padfX[j2] >= cxmin && psCShape->padfX[j2] <= cxmax && + psCShape->padfY[j2] >= cymin && psCShape->padfY[j2] <= cymax ; - if (ierase) inside=(! inside); - if (inside) - { - if (i2 != j2) - { - if (prev_outside) - { - /*** AddIntersection(i2); /*** Add intersection ***/ - prev_outside=FALSE; - } - psCShape->padfX[i2]=psCShape->padfX[j2]; /** move vertex **/ - psCShape->padfY[i2]=psCShape->padfY[j2]; - } - i2++; - } else { - if ( (! prev_outside) && (j2 > 0) ) - { - /*** AddIntersection(i2); /*** Add intersection (Watch out for j2==i2-1) ***/ - /*** Also a polygon may overlap twice and will split into a several parts ***/ - prev_outside=TRUE; - } - } - } + if (ierase) { + inside=(! inside); + } + + if (inside) { + if (i2 != j2) { + if (prev_outside) { + /*** Add intersection ***/ + prev_outside=FALSE; + } + psCShape->padfX[i2]=psCShape->padfX[j2]; /** move vertex **/ + psCShape->padfY[i2]=psCShape->padfY[j2]; + } + i2++; + } else { + if ( (! prev_outside) && (j2 > 0) ) { + /*** Add intersection (Watch out for j2==i2-1) ***/ + /*** Also a polygon may overlap twice and will split into a several parts ***/ + prev_outside=TRUE; + } + } + } printf("Vertices:%d OUT:%d Number of Parts:%d\n", psCShape->nVertices,i2, psCShape->nParts ); - psCShape->nVertices = i2; + psCShape->nVertices = i2; - if (i2 < 2) return(0); /** SKIP RECORD **/ - /*** (WE ARE NOT CREATING INTERESECTIONS and some lines could be reduced to one point) **/ + if (i2 < 2) { + return(0); /** SKIP RECORD **/ + } - if (i2 == 0) return(0); /** SKIP RECORD **/ - else return(1); /** WRITE RECORD **/ - } /** End CUT **/ + /*** (WE ARE NOT CREATING INTERESECTIONS and some lines could be reduced to one point) **/ + if (i2 == 0) { + return(0); /** SKIP RECORD **/ + } else { + return(1); /** WRITE RECORD **/ + } + } /** End CUT **/ + + return -1; } @@ -909,31 +964,29 @@ int j,i; return(0); } - #define NKEYS (sizeof(unitkeytab) / sizeof(struct unitkey)) -findunit(unit) - char *unit; - { +int findunit(char *unit) +{ struct unitkey { char *name; double value; } unitkeytab[] = { - "CM", 39.37, - "CENTIMETER", 39.37, - "CENTIMETERS", 39.37, /** # of inches * 100 in unit **/ - "METER", 3937, - "METERS", 3937, - "KM", 3937000, - "KILOMETER", 3937000, - "KILOMETERS", 3937000, - "INCH", 100, - "INCHES", 100, - "FEET", 1200, - "FOOT", 1200, - "YARD", 3600, - "YARDS", 3600, - "MILE", 6336000, - "MILES", 6336000 + { "CM", 39.37 }, + { "CENTIMETER", 39.37 }, + { "CENTIMETERS", 39.37 }, /** # of inches * 100 in unit **/ + { "METER", 3937 }, + { "METERS", 3937 }, + { "KM", 3937000 }, + { "KILOMETER", 3937000 }, + { "KILOMETERS", 3937000 }, + { "INCH", 100 }, + { "INCHES", 100 }, + { "FEET", 1200 }, + { "FOOT", 1200 }, + { "YARD", 3600 }, + { "YARDS", 3600 }, + { "MILE", 6336000 }, + { "MILES", 6336000 } }; double unitfactor=0; diff --git a/src/Lib/vpf/table.cxx b/src/Lib/vpf/table.cxx index 66dc0023..e3df2425 100644 --- a/src/Lib/vpf/table.cxx +++ b/src/Lib/vpf/table.cxx @@ -529,49 +529,57 @@ VpfTable::read_double (istream &input) const short VpfTable::make_short (char buf[2]) const { - if (_system_byte_order == _file_byte_order) { - return *((short *)buf); - } else { - char out[2]; - swap2(buf, out); - return *((short *)out); - } + if (_system_byte_order == _file_byte_order) { + return *((short *)buf); + } else { + char out[2]; + short *out_short = (short *)out; + + swap2(buf, out); + return *out_short; + } } int VpfTable::make_int (char buf[4]) const { - if (_system_byte_order == _file_byte_order) { - return *((int *)buf); - } else { - char out[4]; - swap4(buf, out); - return *((int *)out); - } + if (_system_byte_order == _file_byte_order) { + return *((int *)buf); + } else { + char out[4]; + int *int_out = (int *)out; + + swap4(buf, out); + return *int_out; + } } float VpfTable::make_float (char buf[4]) const { - if (_system_byte_order == _file_byte_order) { - return *((float *)buf); - } else { - char out[4]; - swap4(buf, out); - return *((float *)out); - } + if (_system_byte_order == _file_byte_order) { + return *((float *)buf); + } else { + char out[4]; + float *float_out = (float *)out; + + swap4(buf, out); + return *float_out; + } } double VpfTable::make_double (char buf[8]) const { - if (_system_byte_order == _file_byte_order) { - return *((double *)buf); - } else { - char out[8]; - swap8(buf, out); - return *((double *)out); - } + if (_system_byte_order == _file_byte_order) { + return *((double *)buf); + } else { + char out[8]; + double *double_out = (double *)out; + + swap8(buf, out); + return *double_out; + } } diff --git a/src/Lib/vpf/value.cxx b/src/Lib/vpf/value.cxx index 4c18e02f..17cee671 100644 --- a/src/Lib/vpf/value.cxx +++ b/src/Lib/vpf/value.cxx @@ -485,37 +485,41 @@ VpfValue::convertType (char rawType) ostream & operator<< (ostream &output, const VpfValue &value) { - switch (value.getType()) { - case VpfValue::TEXT: - output << value.getText(); - break; - case VpfValue::INT: - output << value.getInt(); - break; - case VpfValue::REAL: - output << value.getReal(); - break; - case VpfValue::POINTS: { - int nPoints = value.getElementCount(); - for (int i = 0; i < nPoints; i++) { - if (i > 0) - output << ','; - output << '{' << value.getPoint(i).x << ',' - << value.getPoint(i).y << ',' << value.getPoint(i).z << '}'; + switch (value.getType()) { + case VpfValue::TEXT: + output << value.getText(); + break; + case VpfValue::INT: + output << value.getInt(); + break; + case VpfValue::REAL: + output << value.getReal(); + break; + case VpfValue::POINTS: { + int nPoints = value.getElementCount(); + for (int i = 0; i < nPoints; i++) { + if (i > 0) { + output << ','; + } + output << '{' << value.getPoint(i).x << ',' + << value.getPoint(i).y << ',' + << value.getPoint(i).z << '}'; + } + break; + } + case VpfValue::DATE: + output << value.getDate(); // FIXME + break; + case VpfValue::CROSSREF: + output << value.getCrossRef().current_tile_key << ',' + << value.getCrossRef().next_tile_id << ',' + << value.getCrossRef().next_tile_key; + break; + default: + break; } - break; - } - case VpfValue::DATE: - output << value.getDate(); // FIXME - break; - case VpfValue::CROSSREF: - output << value.getCrossRef().current_tile_key << ',' - << value.getCrossRef().next_tile_id << ',' - << value.getCrossRef().next_tile_key; - break; - } - return output; + return output; } // end of value.cxx diff --git a/src/Prep/GDALChop/gdalchop.cxx b/src/Prep/GDALChop/gdalchop.cxx index c7172f26..e4f81180 100644 --- a/src/Prep/GDALChop/gdalchop.cxx +++ b/src/Prep/GDALChop/gdalchop.cxx @@ -79,33 +79,32 @@ struct SimpleRasterTransformerInfo { double col_step, row_step; }; -int SimpleRasterTransformer(void *pTransformerArg, - int bDstToSrc, int nPointCount, +int SimpleRasterTransformer(void *pTransformerArg, + int bDstToSrc, int nPointCount, double *x, double *y, double *z, int *panSuccess ) { - SimpleRasterTransformerInfo* info=(SimpleRasterTransformerInfo*)pTransformerArg; + SimpleRasterTransformerInfo* info = (SimpleRasterTransformerInfo*)pTransformerArg; int success; - + if (bDstToSrc) { /* transform destination -> source */ - for (int i=0;i<nPointCount;i++) { - x[i]=info->x0+info->col_step*x[i]; - y[i]=info->y0+info->row_step*y[i]; + for (int i = 0; i < nPointCount; i++) { + x[i] = info->x0 + info->col_step * x[i]; + y[i] = info->y0 + info->row_step * y[i]; } - - success=info->pfnTransformer(info->pTransformerArg, - bDstToSrc,nPointCount, - x,y,z,panSuccess); + + success = info->pfnTransformer(info->pTransformerArg, + bDstToSrc, nPointCount, + x, y, z, panSuccess); } else { - success=info->pfnTransformer(info->pTransformerArg, - bDstToSrc,nPointCount, - x,y,z,panSuccess); - for (int i=0;i<nPointCount;i++) { - if (!panSuccess[i]) { + success = info->pfnTransformer(info->pTransformerArg, + bDstToSrc, nPointCount, + x, y, z, panSuccess); + for (int i = 0; i < nPointCount; i++) { + if (!panSuccess[i]) continue; - } - x[i]=(x[i]-info->x0)/info->col_step; - y[i]=(y[i]-info->y0)/info->row_step; + x[i] = (x[i] - info->x0) / info->col_step; + y[i] = (y[i] - info->y0) / info->row_step; } } return success; @@ -113,178 +112,183 @@ int SimpleRasterTransformer(void *pTransformerArg, class ImageInfo { public: - ImageInfo(GDALDataset *dataset); - - void GetBounds(double &n, double &s, double &e, double &w) const { - n=north; - s=south; - e=east; - w=west; - } - - const char* GetDescription() const { - return dataset->GetDescription(); - } - - void GetDataChunk(int *buffer, - double x, double y, - double colstep, double rowstep, - int w, int h, - int srcband=1, int nodata=-32768); +ImageInfo(GDALDataset *dataset); + +void GetBounds(double &n, double &s, double &e, double &w) const +{ + n = north; + s = south; + e = east; + w = west; +} + +const char* GetDescription() const +{ + return dataset->GetDescription(); +} + +void GetDataChunk(int *buffer, + double x, double y, + double colstep, double rowstep, + int w, int h, + int srcband = 1, int nodata = -32768); protected: - /* The dataset */ - GDALDataset *dataset; - - /* Source spatial reference system */ - OGRSpatialReference srs; - - /* Coordinate transformation pixel -> geographic */ - double geoXfrm[6]; - - /* Coordinate transformation to WGS84 */ - OGRCoordinateTransformation *wgs84xform; - - /* geographical edge coordinates in CCW order, WGS84 */ - double geoX[4],geoY[4]; - - /* bounding box in WGS84 */ - double north,south,east,west; +/* The dataset */ +GDALDataset *dataset; + +/* Source spatial reference system */ +OGRSpatialReference srs; + +/* Coordinate transformation pixel -> geographic */ +double geoXfrm[6]; + +/* Coordinate transformation to WGS84 */ +OGRCoordinateTransformation *wgs84xform; + +/* geographical edge coordinates in CCW order, WGS84 */ +double geoX[4], geoY[4]; + +/* bounding box in WGS84 */ +double north, south, east, west; }; -ImageInfo::ImageInfo(GDALDataset *dataset): - dataset(dataset), - srs(dataset->GetProjectionRef()) +ImageInfo::ImageInfo(GDALDataset *dataset) : + dataset(dataset), + srs(dataset->GetProjectionRef()) { OGRSpatialReference wgs84SRS; + wgs84SRS.SetWellKnownGeogCS( "EPSG:4326" ); - - + + /* Determine the bounds of the input file in WGS84 */ - int w=dataset->GetRasterXSize(); - int h=dataset->GetRasterYSize(); - - if (dataset->GetGeoTransform(geoXfrm)!=CE_None) { - SG_LOG(SG_GENERAL, SG_ALERT, - "Could not determine transform matrix for dataset " - "'" << dataset->GetDescription() << "'" - ":" << CPLGetLastErrorMsg()); - exit(1); + int w = dataset->GetRasterXSize(); + int h = dataset->GetRasterYSize(); + + if (dataset->GetGeoTransform(geoXfrm) != CE_None) { + SG_LOG(SG_GENERAL, SG_ALERT, + "Could not determine transform matrix for dataset " + "'" << dataset->GetDescription() << "'" + ":" << CPLGetLastErrorMsg()); + exit(1); } - + /* create points in CCW order */ - geoX[0]=geoXfrm[0]+0*geoXfrm[1]+0*geoXfrm[2]; - geoX[1]=geoXfrm[0]+0*geoXfrm[1]+h*geoXfrm[2]; - geoX[2]=geoXfrm[0]+w*geoXfrm[1]+h*geoXfrm[2]; - geoX[3]=geoXfrm[0]+w*geoXfrm[1]+0*geoXfrm[2]; - geoY[0]=geoXfrm[3]+0*geoXfrm[4]+0*geoXfrm[5]; - geoY[1]=geoXfrm[3]+0*geoXfrm[4]+h*geoXfrm[5]; - geoY[2]=geoXfrm[3]+w*geoXfrm[4]+h*geoXfrm[5]; - geoY[3]=geoXfrm[3]+w*geoXfrm[4]+0*geoXfrm[5]; - - const char* projRef=dataset->GetProjectionRef(); - - srs=OGRSpatialReference(projRef); - + geoX[0] = geoXfrm[0] + 0 * geoXfrm[1] + 0 * geoXfrm[2]; + geoX[1] = geoXfrm[0] + 0 * geoXfrm[1] + h * geoXfrm[2]; + geoX[2] = geoXfrm[0] + w * geoXfrm[1] + h * geoXfrm[2]; + geoX[3] = geoXfrm[0] + w * geoXfrm[1] + 0 * geoXfrm[2]; + geoY[0] = geoXfrm[3] + 0 * geoXfrm[4] + 0 * geoXfrm[5]; + geoY[1] = geoXfrm[3] + 0 * geoXfrm[4] + h * geoXfrm[5]; + geoY[2] = geoXfrm[3] + w * geoXfrm[4] + h * geoXfrm[5]; + geoY[3] = geoXfrm[3] + w * geoXfrm[4] + 0 * geoXfrm[5]; + + const char* projRef = dataset->GetProjectionRef(); + + srs = OGRSpatialReference(projRef); + wgs84xform = OGRCreateCoordinateTransformation( &srs, &wgs84SRS ); - - if (!wgs84xform->Transform(4,geoX,geoY)) { - SG_LOG(SG_GENERAL, SG_ALERT, - "Could not transform edge points of dataset " - "'" << dataset->GetDescription() << "'" - ":" << CPLGetLastErrorMsg()); - exit(1); + + if (!wgs84xform->Transform(4, geoX, geoY)) { + SG_LOG(SG_GENERAL, SG_ALERT, + "Could not transform edge points of dataset " + "'" << dataset->GetDescription() << "'" + ":" << CPLGetLastErrorMsg()); + exit(1); } - - east=west=geoX[0]; - north=south=geoY[0]; - - for (int j=1;j<4;j++) { - north=std::max(north,geoY[j]); - south=std::min(south,geoY[j]); - east =std::max(east ,geoX[j]); - west =std::min(west ,geoX[j]); + + east = west = geoX[0]; + north = south = geoY[0]; + + for (int j = 1; j < 4; j++) { + north = std::max(north, geoY[j]); + south = std::min(south, geoY[j]); + east = std::max(east, geoX[j]); + west = std::min(west, geoX[j]); } - + SG_LOG(SG_GENERAL, SG_INFO, - "INFO: Bounds for '" << dataset->GetDescription() << "' are" - " n=" << north << " s=" << south << - " e=" << east << " w=" << west); + "INFO: Bounds for '" << dataset->GetDescription() << "' are" + " n=" << north << " s=" << south << + " e=" << east << " w=" << west); } void ImageInfo::GetDataChunk(int *buffer, double x, double y, double colstep, double rowstep, int w, int h, - int srcband, int nodata) { + int srcband, int nodata) +{ OGRSpatialReference wgs84SRS; + wgs84SRS.SetWellKnownGeogCS( "EPSG:4326" ); - + char* wgs84WKT; - + wgs84SRS.exportToWkt(&wgs84WKT); - + /* Setup a raster transformation from WGS84 to raster coordinates of the array files */ SimpleRasterTransformerInfo xformData; - xformData.pTransformerArg=GDALCreateGenImgProjTransformer( - dataset,NULL, - NULL,wgs84WKT, - FALSE, - 0.0, - 1); - xformData.pfnTransformer=GDALGenImgProjTransform; - xformData.x0=x; - xformData.y0=y; - xformData.col_step=colstep; - xformData.row_step=rowstep; - + xformData.pTransformerArg = GDALCreateGenImgProjTransformer( + dataset, NULL, + NULL, wgs84WKT, + FALSE, + 0.0, + 1); + xformData.pfnTransformer = GDALGenImgProjTransform; + xformData.x0 = x; + xformData.y0 = y; + xformData.col_step = colstep; + xformData.row_step = rowstep; + // TODO: check if this image can actually cover part of the chunk - + /* establish the full source to target transformation */ GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); - - int srcBandNumbers[]={srcband}; - int dstBandNumbers[]={1}; - double dstNodata[]={nodata}; + + int srcBandNumbers[] = { srcband }; + int dstBandNumbers[] = { 1 }; + double dstNodata[] = { nodata }; double srcNodataReal[1]; - double srcNodataImag[]={0.0}; + double srcNodataImag[] = { 0.0 }; int srcHasNodataValue; - - srcNodataReal[0]=dataset->GetRasterBand(srcband)->GetNoDataValue(&srcHasNodataValue); - - psWarpOptions->hSrcDS=dataset; - psWarpOptions->hDstDS=NULL; - psWarpOptions->nBandCount=1; - psWarpOptions->panSrcBands=srcBandNumbers; - psWarpOptions->panDstBands=dstBandNumbers; - psWarpOptions->nSrcAlphaBand=0; - psWarpOptions->nDstAlphaBand=0; - psWarpOptions->padfSrcNoDataReal=(srcHasNodataValue?srcNodataReal:NULL); - psWarpOptions->padfSrcNoDataImag=(srcHasNodataValue?srcNodataImag:NULL); - psWarpOptions->padfDstNoDataReal=dstNodata; + + srcNodataReal[0] = dataset->GetRasterBand(srcband)->GetNoDataValue(&srcHasNodataValue); + + psWarpOptions->hSrcDS = dataset; + psWarpOptions->hDstDS = NULL; + psWarpOptions->nBandCount = 1; + psWarpOptions->panSrcBands = srcBandNumbers; + psWarpOptions->panDstBands = dstBandNumbers; + psWarpOptions->nSrcAlphaBand = 0; + psWarpOptions->nDstAlphaBand = 0; + psWarpOptions->padfSrcNoDataReal = (srcHasNodataValue ? srcNodataReal : NULL); + psWarpOptions->padfSrcNoDataImag = (srcHasNodataValue ? srcNodataImag : NULL); + psWarpOptions->padfDstNoDataReal = dstNodata; psWarpOptions->eResampleAlg = GRA_Bilinear; psWarpOptions->eWorkingDataType = GDT_Int32; - - psWarpOptions->pfnTransformer=SimpleRasterTransformer; - psWarpOptions->pTransformerArg=&xformData; - + + psWarpOptions->pfnTransformer = SimpleRasterTransformer; + psWarpOptions->pTransformerArg = &xformData; + GDALWarpOperation oOperation; - + oOperation.Initialize( psWarpOptions ); - + /* do the warp */ - if (oOperation.WarpRegionToBuffer(0,0,w,h,buffer,GDT_Int32)!=CE_None) { + if (oOperation.WarpRegionToBuffer(0, 0, w, h, buffer, GDT_Int32) != CE_None) { SG_LOG(SG_GENERAL, SG_ALERT, - "Could not warp to buffer on dataset '" << GetDescription() << "'" - ":" << CPLGetLastErrorMsg()); + "Could not warp to buffer on dataset '" << GetDescription() << "'" + ":" << CPLGetLastErrorMsg()); } - + /* clean up */ - psWarpOptions->panSrcBands=NULL; - psWarpOptions->panDstBands=NULL; - psWarpOptions->padfSrcNoDataReal=NULL; - psWarpOptions->padfSrcNoDataImag=NULL; - psWarpOptions->padfDstNoDataReal=NULL; - + psWarpOptions->panSrcBands = NULL; + psWarpOptions->panDstBands = NULL; + psWarpOptions->padfSrcNoDataReal = NULL; + psWarpOptions->padfSrcNoDataImag = NULL; + psWarpOptions->padfDstNoDataReal = NULL; + GDALDestroyGenImgProjTransformer( xformData.pTransformerArg ); GDALDestroyWarpOptions( psWarpOptions ); } @@ -294,117 +298,117 @@ void write_bucket(const string& work_dir, SGBucket bucket, int min_x, int min_y, int span_x, int span_y, int col_step, int row_step, - bool compress=true) { - SGPath path(work_dir); - path.append(bucket.gen_base_path()); - path.create_dir( 0755 ); + bool compress = true) +{ + SGPath path(work_dir); - string array_file = path.str() + "/" + bucket.gen_index_str() + ".arr"; - - FILE *fp; - if ( (fp = fopen(array_file.c_str(), "w")) == NULL ) { - SG_LOG(SG_GENERAL, SG_ALERT, "cannot open " << array_file << " for writing!"); - exit(-1); - } - - fprintf( fp, "%d %d\n", min_x, min_y ); - fprintf( fp, "%d %d %d %d\n", - span_x + 1, col_step, - span_y + 1, row_step ); - int k=0; - for ( int x = 0; x <= span_x; ++x ) { - for ( int y = 0; y <= span_y; ++y ) { - fprintf( fp, "%d ", buffer[ y * span_x + x ] ); - } - fprintf( fp, "\n" ); - } - fclose(fp); - - if ( compress ) { - string command = "gzip --best -f " + array_file; - system( command.c_str() ); - } + path.append(bucket.gen_base_path()); + path.create_dir( 0755 ); + + string array_file = path.str() + "/" + bucket.gen_index_str() + ".arr"; + + FILE *fp; + if ( (fp = fopen(array_file.c_str(), "w")) == NULL ) { + SG_LOG(SG_GENERAL, SG_ALERT, "cannot open " << array_file << " for writing!"); + exit(-1); + } + + fprintf( fp, "%d %d\n", min_x, min_y ); + fprintf( fp, "%d %d %d %d\n", + span_x + 1, col_step, + span_y + 1, row_step ); + + for ( int x = 0; x <= span_x; ++x ) { + for ( int y = 0; y <= span_y; ++y ) + fprintf( fp, "%d ", buffer[ y * span_x + x ] ); + fprintf( fp, "\n" ); + } + fclose(fp); + + if ( compress ) { + string command = "gzip --best -f " + array_file; + system( command.c_str() ); + } } void process_bucket(const string& work_dir, SGBucket bucket, ImageInfo* images[], int imagecount, - bool forceWrite=false) { - double clat,clon; - double bnorth,bsouth,beast,bwest; - - clat=bucket.get_center_lat(); - clon=bucket.get_center_lon(); - - bnorth=clat+bucket.get_height()/2.0; - bsouth=clat-bucket.get_height()/2.0; - beast =clon+bucket.get_width()/2.0; - bwest =clon-bucket.get_width()/2.0; - - SG_LOG(SG_GENERAL, SG_INFO, "processing bucket " << bucket << "(" << bucket.gen_index() << ") n=" << bnorth << " s=" << bsouth << " e=" << beast << " w=" << bwest); - - /* Get the data from the input datasets... */ - int min_x, min_y, span_x, span_y; - - min_x = (int)(bwest*3600.0); - min_y = (int)(bsouth*3600.0); - - // TODO: Make other resolutions possible as well - int col_step=3, row_step=3; - span_x=bucket.get_width()*3600/col_step; - span_y=bucket.get_height()*3600/row_step; - - int cellcount=(span_x+1)*(span_y+1); - boost::scoped_array<int> buffer(new int[cellcount]); - - ::memset(buffer.get(),-1,(span_x+1)*(span_y+1)*sizeof(int)); - - for (int i=0;i<imagecount;i++) { - double inorth,isouth,ieast,iwest; - images[i]->GetBounds(inorth,isouth,ieast,iwest); - - images[i]->GetDataChunk(buffer.get(), - bwest, bsouth, - col_step/3600.0, row_step/3600.0, - span_x+1, span_y+1); - } - - /* ...check the amount of undefined cells... */ - int nodataCellCount=0; - - for (int i=0;i<cellcount;i++) { - if (buffer[i]==-32768) { - nodataCellCount++; - } - } - - const double nodataPercLimit=5.0; - double nodataPerc=100.0*nodataCellCount/cellcount; - - if (nodataCellCount>0) { - SG_LOG(SG_GENERAL, SG_INFO, " " << nodataCellCount << " cells are not covered with data (" << nodataPerc << "% of cells)"); - } - if (nodataPerc>nodataPercLimit) { - SG_LOG(SG_GENERAL, SG_INFO, " there is not enough data available to cover this cell (limit for non-covered cells is " << nodataPercLimit <<"%)"); - /* don't write out if not forced to */ - if (!forceWrite) - return; - } - - /* ...and write it out */ - write_bucket(work_dir, bucket, - buffer.get(), - min_x, min_y, - span_x, span_y, - col_step, row_step); + bool forceWrite = false) +{ + double clat, clon; + double bnorth, bsouth, beast, bwest; + + clat = bucket.get_center_lat(); + clon = bucket.get_center_lon(); + + bnorth = clat + bucket.get_height() / 2.0; + bsouth = clat - bucket.get_height() / 2.0; + beast = clon + bucket.get_width() / 2.0; + bwest = clon - bucket.get_width() / 2.0; + + SG_LOG(SG_GENERAL, SG_INFO, "processing bucket " << bucket << "(" << bucket.gen_index() << ") n=" << bnorth << " s=" << bsouth << " e=" << beast << " w=" << bwest); + + /* Get the data from the input datasets... */ + int min_x, min_y, span_x, span_y; + + min_x = (int)(bwest * 3600.0); + min_y = (int)(bsouth * 3600.0); + + // TODO: Make other resolutions possible as well + int col_step = 3, row_step = 3; + span_x = bucket.get_width() * 3600 / col_step; + span_y = bucket.get_height() * 3600 / row_step; + + int cellcount = (span_x + 1) * (span_y + 1); + boost::scoped_array<int> buffer(new int[cellcount]); + + ::memset(buffer.get(), -1, (span_x + 1) * (span_y + 1) * sizeof(int)); + + for (int i = 0; i < imagecount; i++) { + double inorth, isouth, ieast, iwest; + images[i]->GetBounds(inorth, isouth, ieast, iwest); + + images[i]->GetDataChunk(buffer.get(), + bwest, bsouth, + col_step / 3600.0, row_step / 3600.0, + span_x + 1, span_y + 1); + } + + /* ...check the amount of undefined cells... */ + int nodataCellCount = 0; + + for (int i = 0; i < cellcount; i++) + if (buffer[i] == -32768) + nodataCellCount++; + + const double nodataPercLimit = 5.0; + double nodataPerc = 100.0 * nodataCellCount / cellcount; + + if (nodataCellCount > 0) + SG_LOG(SG_GENERAL, SG_INFO, " " << nodataCellCount << " cells are not covered with data (" << nodataPerc << "% of cells)"); + if (nodataPerc > nodataPercLimit) { + SG_LOG(SG_GENERAL, SG_INFO, " there is not enough data available to cover this cell (limit for non-covered cells is " << nodataPercLimit << "%)"); + /* don't write out if not forced to */ + if (!forceWrite) + return; + } + + /* ...and write it out */ + write_bucket(work_dir, bucket, + buffer.get(), + min_x, min_y, + span_x, span_y, + col_step, row_step); } -int main(int argc, const char **argv) { +int main(int argc, const char **argv) +{ sglog().setLogLevels( SG_ALL, SG_INFO ); if ( argc < 3 ) { - SG_LOG(SG_GENERAL, SG_ALERT, - "Usage " << argv[0] << " <work_dir> <datasetname...> [-- <bucket-idx> ...]"); - exit(-1); + SG_LOG(SG_GENERAL, SG_ALERT, + "Usage " << argv[0] << " <work_dir> <datasetname...> [-- <bucket-idx> ...]"); + exit(-1); } SGPath work_dir( argv[1] ); @@ -413,61 +417,59 @@ int main(int argc, const char **argv) { GDALAllRegister(); - int datasetcount=0,tilecount=0; + int datasetcount = 0, tilecount = 0; int dashpos; - for (dashpos=2;dashpos<argc;dashpos++) { - if (!strcmp(argv[dashpos],"--")) { - break; - } + for (dashpos = 2; dashpos < argc; dashpos++) + if (!strcmp(argv[dashpos], "--")) + break; + + datasetcount = dashpos - 2; + tilecount = (dashpos == argc ? 0 : argc - dashpos - 1); + + if (datasetcount == 0) { + SG_LOG(SG_GENERAL, SG_ALERT, + "No data sets supplied. Must provide at least one dataset."); + exit(1); } - datasetcount=dashpos-2; - tilecount=(dashpos==argc?0:argc-dashpos-1); + const char** tilenames = argv + dashpos + 1; + const char** datasetnames = argv + 2; - if (datasetcount==0) { - SG_LOG(SG_GENERAL, SG_ALERT, - "No data sets supplied. Must provide at least one dataset."); - exit(1); - } - - const char** tilenames=argv+dashpos+1; - const char** datasetnames=argv+2; - boost::scoped_array<ImageInfo *> images( new ImageInfo *[datasetcount] ); - - double north=-1000,south=1000,east=-1000,west=1000; - + + double north = -1000, south = 1000, east = -1000, west = 1000; + // TODO: allow specification of bounds by user - + /* * Step 1: Open all provided datasets and determine their bounds in WGS84. */ - for (int i=0;i<datasetcount;i++) { - GDALDataset* dataset; - - dataset=(GDALDataset *)GDALOpenShared(datasetnames[i],GA_ReadOnly); - - if (dataset==NULL) { - SG_LOG(SG_GENERAL, SG_ALERT, - "Could not open dataset '" << datasetnames[i] << "'" - ":" << CPLGetLastErrorMsg()); - exit(1); - } - - images[i]=new ImageInfo(dataset); - - double inorth,isouth,ieast,iwest; - images[i]->GetBounds(inorth,isouth,ieast,iwest); - - north=std::max(north,inorth); - south=std::min(south,isouth); - east =std::max(east ,ieast ); - west =std::min(west ,iwest ); + for (int i = 0; i < datasetcount; i++) { + GDALDataset* dataset; + + dataset = (GDALDataset*)GDALOpenShared(datasetnames[i], GA_ReadOnly); + + if (dataset == NULL) { + SG_LOG(SG_GENERAL, SG_ALERT, + "Could not open dataset '" << datasetnames[i] << "'" + ":" << CPLGetLastErrorMsg()); + exit(1); + } + + images[i] = new ImageInfo(dataset); + + double inorth, isouth, ieast, iwest; + images[i]->GetBounds(inorth, isouth, ieast, iwest); + + north = std::max(north, inorth); + south = std::min(south, isouth); + east = std::max(east, ieast ); + west = std::min(west, iwest ); } - + SG_LOG(SG_GENERAL, SG_INFO, "Complete bounds n=" << north << " s=" << south << " e=" << east << " w=" << west); - + /* * Step 2: If no tiles were specified, go through all tiles contained in * the common bounds of all datasets and find those which have @@ -476,25 +478,25 @@ int main(int argc, const char **argv) { * all of them. Warn if no sufficient coverage (non-null pixels) is * available. */ - if (tilecount==0) { + if (tilecount == 0) { /* * No tiles were specified, so we determine the common bounds of all * specified datasets and check all the tiles contained in them. */ - - SGBucket start(west,south),end(east,north); - - int dx,dy; - + + SGBucket start(west, south), end(east, north); + + int dx, dy; + sgBucketDiff(start, end, &dx, &dy); - + SG_LOG(SG_GENERAL, SG_INFO, "dx=" << dx << " dy=" << dy); - - for (int x=0;x<=dx;x++) { - for (int y=0;y<=dy;y++) { - SGBucket bucket=sgBucketOffset(west,south,x,y); - - process_bucket(work_dir.str(),bucket,images.get(),datasetcount); + + for (int x = 0; x <= dx; x++) { + for (int y = 0; y <= dy; y++) { + SGBucket bucket = sgBucketOffset(west, south, x, y); + + process_bucket(work_dir.str(), bucket, images.get(), datasetcount); } } } else { @@ -502,13 +504,13 @@ int main(int argc, const char **argv) { * Tiles were specified, so process them and warn if not enough * data is available, but write them in any case. */ - for (int i=0;i<tilecount;i++) { + for (int i = 0; i < tilecount; i++) { SGBucket bucket(atol(tilenames[i])); - - process_bucket(work_dir.str(),bucket,images.get(),datasetcount,true); + + process_bucket(work_dir.str(), bucket, images.get(), datasetcount, true); } } - + return 0; } diff --git a/src/Prep/Terra/Subdivision.cc b/src/Prep/Terra/Subdivision.cc index 12cda725..b34fc358 100644 --- a/src/Prep/Terra/Subdivision.cc +++ b/src/Prep/Terra/Subdivision.cc @@ -139,7 +139,7 @@ static unsigned int timestamp = 0; static void overEdge(Edge *e, edge_callback fn, void *closure) { - if( e->token != timestamp ) + if( e->token != (int)timestamp ) { e->token = timestamp; e->Sym()->token = timestamp; @@ -247,37 +247,42 @@ Edge *Subdivision::locate(const Vec2& x, Edge *start) real to = triArea(x, eo->Dest(), eo->Org()); real td = triArea(x, ed->Dest(), ed->Org()); - if (td>0) // x is below ed - if (to>0 || to==0 && t==0) {// x is interior, or origin endpoint + if ( td > 0 ) { // x is below ed + if ( (to > 0) || (to==0 && t==0) ) { // x is interior, or origin endpoint startingEdge = e; return e; } - else { // x is below ed, below eo + else { // x is below ed, below eo t = to; e = eo; } - else // x is on or above ed - if (to>0) // x is above eo - if (td==0 && t==0) { // x is destination endpoint + } + else { // x is on or above ed + if (to>0) { // x is above eo + if (td==0 && t==0) { // x is destination endpoint startingEdge = e; return e; } - else { // x is on or above ed and above eo + else { // x is on or above ed and above eo t = td; e = ed; } - else // x is on or below eo - if (t==0 && !leftOf(eo->Dest(), e)) + } + else { // x is on or below eo + if (t==0 && !leftOf(eo->Dest(), e)) { // x on e but subdiv. is to right e = e->Sym(); - else if (rand()&1) { // x is on or above ed and - t = to; // on or below eo; step randomly + } + else if (rand()&1) { // x is on or above ed and + t = to; // on or below eo; step randomly e = eo; } else { t = td; e = ed; } + } + } } } @@ -372,7 +377,7 @@ void Subdivision::optimize(Vec2& x, Edge *s) do { Edge *e = spoke->Lnext(); - Edge *t = e->Oprev(); + // Edge *t = e->Oprev(); if( isInterior(e) && shouldSwap(x, e) ) swap(e);