diff --git a/src/Lib/terragear/clipper.cpp b/src/Lib/terragear/clipper.cpp index e3e61786..060ab6c8 100644 --- a/src/Lib/terragear/clipper.cpp +++ b/src/Lib/terragear/clipper.cpp @@ -1,8 +1,8 @@ /******************************************************************************* * * * Author : Angus Johnson * -* Version : 5.1.2 * -* Date : 25 February 2013 * +* Version : 5.1.4 * +* Date : 24 March 2013 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2013 * * * @@ -543,14 +543,18 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2, IntPoint &ip, bool UseFullInt64Range) { double b1, b2; - if (SlopesEqual(edge1, edge2, UseFullInt64Range)) return false; + if (SlopesEqual(edge1, edge2, UseFullInt64Range)) + { + if (edge2.ybot > edge1.ybot) ip.Y = edge2.ybot; + else ip.Y = edge1.ybot; + return false; + } else if (NEAR_ZERO(edge1.dx)) { ip.X = edge1.xbot; if (NEAR_EQUAL(edge2.dx, HORIZONTAL)) - { ip.Y = edge2.ybot; - } else + else { b2 = edge2.ybot - (edge2.xbot / edge2.dx); ip.Y = Round(ip.X / edge2.dx + b2); @@ -560,14 +564,14 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2, { ip.X = edge2.xbot; if (NEAR_EQUAL(edge1.dx, HORIZONTAL)) - { ip.Y = edge1.ybot; - } else + else { b1 = edge1.ybot - (edge1.xbot / edge1.dx); ip.Y = Round(ip.X / edge1.dx + b1); } - } else + } + else { b1 = edge1.xbot - edge1.ybot * edge1.dx; b2 = edge2.xbot - edge2.ybot * edge2.dx; @@ -586,7 +590,8 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2, ip.X = edge1.xtop; ip.Y = edge1.ytop; return TopX(edge2, edge1.ytop) < edge1.xtop; - } else + } + else { ip.X = edge2.xtop; ip.Y = edge2.ytop; @@ -1761,17 +1766,11 @@ void Clipper::IntersectEdges(TEdge *e1, TEdge *e2, } else if ( e1Contributing ) { - if ((e2Wc == 0 || e2Wc == 1) && - (m_ClipType != ctIntersection || - e2->polyType == ptSubject || (e2->windCnt2 != 0))) - DoEdge1(e1, e2, pt); + if (e2Wc == 0 || e2Wc == 1) DoEdge1(e1, e2, pt); } else if ( e2contributing ) { - if ((e1Wc == 0 || e1Wc == 1) && - (m_ClipType != ctIntersection || - e1->polyType == ptSubject || (e1->windCnt2 != 0))) - DoEdge2(e1, e2, pt); + if (e1Wc == 0 || e1Wc == 1) DoEdge2(e1, e2, pt); } else if ( (e1Wc == 0 || e1Wc == 1) && (e2Wc == 0 || e2Wc == 1) && !e1stops && !e2stops ) @@ -2204,28 +2203,28 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge) TEdge* e = GetNextInAEL( horzEdge , dir ); while( e ) { + if ( e->xcurr == horzEdge->xtop && !eMaxPair ) + { + if (SlopesEqual(*e, *horzEdge->nextInLML, m_UseFullRange)) + { + //if output polygons share an edge, they'll need joining later ... + if (horzEdge->outIdx >= 0 && e->outIdx >= 0) + AddJoin(horzEdge->nextInLML, e, horzEdge->outIdx); + break; //we've reached the end of the horizontal line + } + else if (e->dx < horzEdge->nextInLML->dx) + //we really have got to the end of the intermediate horz edge so quit. + //nb: More -ve slopes follow more +ve slopes ABOVE the horizontal. + break; + } + TEdge* eNext = GetNextInAEL( e, dir ); if (eMaxPair || - ((dir == dLeftToRight) && (e->xcurr <= horzRight)) || - ((dir == dRightToLeft) && (e->xcurr >= horzLeft))) + ((dir == dLeftToRight) && (e->xcurr < horzRight)) || + ((dir == dRightToLeft) && (e->xcurr > horzLeft))) { - //ok, so far it looks like we're still in range of the horizontal edge - if ( e->xcurr == horzEdge->xtop && !eMaxPair ) - { - if (SlopesEqual(*e, *horzEdge->nextInLML, m_UseFullRange)) - { - //if output polygons share an edge, they'll need joining later ... - if (horzEdge->outIdx >= 0 && e->outIdx >= 0) - AddJoin(horzEdge->nextInLML, e, horzEdge->outIdx); - break; //we've reached the end of the horizontal line - } - else if (e->dx < horzEdge->nextInLML->dx) - //we really have got to the end of the intermediate horz edge so quit. - //nb: More -ve slopes follow more +ve slopes ABOVE the horizontal. - break; - } - + //so far we're still in range of the horizontal edge if( e == eMaxPair ) { //horzEdge is evidently a maxima horizontal and we've arrived at its end. @@ -2261,8 +2260,8 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge) } SwapPositionsInAEL( horzEdge, e ); } - else if( (dir == dLeftToRight && e->xcurr > horzRight && m_SortedEdges) || - (dir == dRightToLeft && e->xcurr < horzLeft && m_SortedEdges) ) break; + else if( (dir == dLeftToRight && e->xcurr >= horzRight) || + (dir == dRightToLeft && e->xcurr <= horzLeft) ) break; e = eNext; } //end while @@ -2345,7 +2344,7 @@ void Clipper::BuildIntersectList(const long64 botY, const long64 topY) { e->prevInSEL = e->prevInAEL; e->nextInSEL = e->nextInAEL; - e->tmpX = TopX( *e, topY ); + e->xcurr = TopX( *e, topY ); e = e->nextInAEL; } @@ -2359,9 +2358,10 @@ void Clipper::BuildIntersectList(const long64 botY, const long64 topY) { TEdge *eNext = e->nextInSEL; IntPoint pt; - if(e->tmpX > eNext->tmpX && - IntersectPoint(*e, *eNext, pt, m_UseFullRange)) + if(e->xcurr > eNext->xcurr) { + if (!IntersectPoint(*e, *eNext, pt, m_UseFullRange) && e->xcurr > eNext->xcurr +1) + throw clipperException("Intersection error"); if (pt.Y > botY) { pt.Y = botY; @@ -2452,7 +2452,7 @@ void Clipper::DoMaxima(TEdge *e, long64 topY) if (!eNext) throw clipperException("DoMaxima error"); IntersectEdges( e, eNext, IntPoint(X, topY), ipBoth ); SwapPositionsInAEL(e, eNext); - eNext = eNext->nextInAEL; + eNext = e->nextInAEL; } if( e->outIdx < 0 && eMaxPair->outIdx < 0 ) { @@ -2600,25 +2600,22 @@ void Clipper::FixupOutPolygon(OutRec &outRec) void Clipper::BuildResult(Polygons &polys) { - int k = 0; - polys.resize(m_PolyOuts.size()); + polys.reserve(m_PolyOuts.size()); for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); ++i) { if (m_PolyOuts[i]->pts) { - Polygon* pg = &polys[k]; - pg->clear(); + Polygon pg; OutPt* p = m_PolyOuts[i]->pts; do { - pg->push_back(p->pt); + pg.push_back(p->pt); p = p->prev; } while (p != m_PolyOuts[i]->pts); - //make sure each polygon has at least 3 vertices ... - if (pg->size() < 3) pg->clear(); else k++; + if (pg.size() > 2) + polys.push_back(pg); } } - polys.resize(k); } //------------------------------------------------------------------------------ @@ -3178,14 +3175,14 @@ PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys, switch (jointype) { - case jtRound: //limit defaults to 0.125 - if (limit <= 0) limit = 0.125; + case jtRound: + if (limit <= 0) limit = 0.25; else if (limit > std::fabs(delta)) limit = std::fabs(delta); break; - case jtMiter: //limit defaults to twice delta's width ... + case jtMiter: if (limit < 2) limit = 2; break; - default: //otherwise limit is unused + default: //unused limit = 1; } m_RMin = 2.0/(limit*limit); @@ -3275,9 +3272,8 @@ private: void AddPoint(const IntPoint& pt) { - Polygon::size_type len = m_curr_poly->size(); - if (len == m_curr_poly->capacity()) - m_curr_poly->reserve(len + buffLength); + if (m_curr_poly->size() == m_curr_poly->capacity()) + m_curr_poly->reserve(m_curr_poly->capacity() + buffLength); m_curr_poly->push_back(pt); } //------------------------------------------------------------------------------ @@ -3410,38 +3406,62 @@ void SimplifyPolygons(Polygons &polys, PolyFillType fillType) } //------------------------------------------------------------------------------ -void CleanPolygon(Polygon& in_poly, Polygon& out_poly, double distance) +bool PointsAreClose(IntPoint pt1, IntPoint pt2, long64 distSqrd) { - //delta = proximity in units/pixels below which vertices - //will be stripped. Default ~= sqrt(2) so when adjacent - //vertices have both x & y coords within 1 unit, then - //the second vertex will be stripped. - int len = in_poly.size(); - if (len < 3) - out_poly.resize(0); - else - out_poly.resize(in_poly.size()); - - int d = (int)(distance * distance); - IntPoint p = in_poly[0]; - int j = 1; - for (int i = 1; i < len; i++) - { - if ((in_poly[i].X - p.X) * (in_poly[i].X - p.X) + - (in_poly[i].Y - p.Y) * (in_poly[i].Y - p.Y) <= d) - continue; - out_poly[j] = in_poly[i]; - p = in_poly[i]; - j++; - } - p = in_poly[j - 1]; - if ((in_poly[0].X - p.X) * (in_poly[0].X - p.X) + - (in_poly[0].Y - p.Y) * (in_poly[0].Y - p.Y) <= d) - j--; - if (j < len) - out_poly.resize(j); + long64 dx = pt1.X - pt2.X; + long64 dy = pt1.Y - pt2.Y; + return ((dx * dx) + (dy * dy) <= distSqrd); } //------------------------------------------------------------------------------ + +void CleanPolygon(Polygon& in_poly, Polygon& out_poly, double distance) +{ + //distance = proximity in units/pixels below which vertices + //will be stripped. Default ~= sqrt(2). + int highI = in_poly.size() -1; + long64 d = (int)(distance * distance); + while (highI > 0 && PointsAreClose(in_poly[highI], in_poly[0], d)) highI--; + if (highI < 2) + { + out_poly.clear(); + return; + } + out_poly.resize(highI + 1); + bool UseFullRange = FullRangeNeeded(in_poly); + IntPoint pt = in_poly[highI]; + int i = 0; + int k = 0; + for (;;) + { + if (i >= highI) break; + int j = i + 1; + + if (PointsAreClose(pt, in_poly[j], d)) + { + i = j + 1; + while (i <= highI && PointsAreClose(pt, in_poly[i], d)) i++; + continue; + } + + if (PointsAreClose(in_poly[i], in_poly[j], d) || + SlopesEqual(pt, in_poly[i], in_poly[j], UseFullRange)) + { + i = j; + continue; + } + + pt = in_poly[i++]; + out_poly[k++] = pt; + } + + if (i <= highI) out_poly[k++] = in_poly[i]; + if (k > 2 && SlopesEqual(out_poly[k -2], out_poly[k -1], out_poly[0], UseFullRange)) + k--; + if (k < 3) out_poly.clear(); + else if (k <= highI) out_poly.resize(k); +} +//------------------------------------------------------------------------------ + void CleanPolygons(Polygons& in_polys, Polygons& out_polys, double distance) { for (Polygons::size_type i = 0; i < in_polys.size(); ++i) diff --git a/src/Lib/terragear/clipper.hpp b/src/Lib/terragear/clipper.hpp index 1910ec7f..16540f32 100644 --- a/src/Lib/terragear/clipper.hpp +++ b/src/Lib/terragear/clipper.hpp @@ -1,8 +1,8 @@ /******************************************************************************* * * * Author : Angus Johnson * -* Version : 5.1.2 * -* Date : 25 February 2013 * +* Version : 5.1.4 * +* Date : 24 March 2013 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2013 * * * @@ -134,7 +134,6 @@ struct TEdge { double dx; long64 deltaX; long64 deltaY; - long64 tmpX; PolyType polyType; EdgeSide side; int windDelta; //1 or -1 depending on winding direction diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx index fba63ee1..092678bf 100644 --- a/src/Lib/terragear/tg_chopper.cxx +++ b/src/Lib/terragear/tg_chopper.cxx @@ -8,6 +8,9 @@ #include "tg_chopper.hxx" +#include "tg_shapefile.hxx" + +unsigned int strip_clip = 0; tgPolygon tgChopper::Clip( const tgPolygon& subject, const std::string& type, @@ -42,7 +45,7 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject, base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), max.getLatitudeDeg()) ); base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), max.getLatitudeDeg()) ); - result = tgPolygon::Intersect( subject, base ); + result = tgPolygon::Intersect( base, subject ); if ( result.Contours() > 0 ) { if ( subject.GetPreserve3D() ) { result.InheritElevations( subject ); @@ -63,17 +66,33 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject, return result; } +// Pass in the center lat for clipping buckets from the row. +// We can't rely on sgBucketOffset, as rounding error sometimes causes it to look like there are 2 rows +// (the first being a sliver) +// This leads to using that poly as the subject - which leads to having no usable polygon for this row. +void tgChopper::ClipRow( const tgPolygon& subject, const double& center_lat, const std::string& type ) +{ + tgRectangle bb = subject.GetBoundingBox(); + SGBucket b_min( bb.getMin() ); + SGBucket b_max( bb.getMax() ); + double min_center_lon = b_min.get_center_lon(); + int dx, dy; + + sgBucketDiff(b_min, b_max, &dx, &dy); + + for ( int i = 0; i <= dx; ++i ) { + SGBucket b_cur = sgBucketOffset(min_center_lon, center_lat, i, 0); + Clip( subject, type, b_cur ); + } +} + void tgChopper::Add( const tgPolygon& subject, const std::string& type ) { - tgRectangle bb; - SGGeod p; - // bail out immediately if polygon is empty if ( subject.Contours() == 0 ) return; - bb = subject.GetBoundingBox(); - + tgRectangle bb = subject.GetBoundingBox(); SG_LOG( SG_GENERAL, SG_DEBUG, " min = " << bb.getMin() << " max = " << bb.getMax() ); // find buckets for min, and max points of convex hull. @@ -81,99 +100,58 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type ) // polygons that span the date line SGBucket b_min( bb.getMin() ); SGBucket b_max( bb.getMax() ); - SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket min = " << b_min ); - SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket max = " << b_max ); + SGBucket b_cur; + int dx, dy; - if ( b_min == b_max ) { - // shape entirely contained in a single bucket, write and bail - Clip( subject, type, b_min ); - } else { - SGBucket b_cur; - int dx, dy; + sgBucketDiff(b_min, b_max, &dx, &dy); + SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy ); - sgBucketDiff(b_min, b_max, &dx, &dy); - SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" ); - SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy ); + if ( (dx > 2880) || (dy > 1440) ) + throw sg_exception("something is really wrong in split_polygon()!!!!"); - if ( (dx > 2880) || (dy > 1440) ) - throw sg_exception("something is really wrong in split_polygon()!!!!"); + if ( dy == 0 ) + { + // We just have a single row - no need to intersect first + ClipRow( subject, b_min.get_center_lat(), type ); + } + else + { + // Multiple rows - perform row intersection to reduce the number of bucket clips we need + // since many shapes are narraw in some places, wide in others - bb will be at the widest part + SG_LOG( SG_GENERAL, SG_DEBUG, "subject spans tile rows: bb is from lat " << bb.getMin().getLatitudeDeg() << " to " << bb.getMax().getLatitudeDeg() << " dy is " << dy ); - if ( dy <= 1 ) { - // we are down to at most two rows, write each column and then bail - double min_center_lat = b_min.get_center_lat(); - double min_center_lon = b_min.get_center_lon(); - for ( int j = 0; j <= dy; ++j ) { - for ( int i = 0; i <= dx; ++i ) { - b_cur = sgBucketOffset(min_center_lon, min_center_lat, i, j); - Clip( subject, type, b_cur ); + for ( int row = 0; row <= dy; row++ ) + { + // Generate a clip rectangle - add some buffer on top and bottom, so we don't clip directly on an edge when we + // clip the individual buckets + // TODO : May no longer be necessary + SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, row ); + double clip_bottom = b_clip.get_center_lat() - SG_HALF_BUCKET_SPAN; // + 0.01); + double clip_top = b_clip.get_center_lat() + SG_HALF_BUCKET_SPAN; // + 0.01); + tgPolygon clip_row, clipped; + + SG_LOG( SG_GENERAL, SG_DEBUG, " row " << row << " center lat is " << b_clip.get_center_lat() << " clip_botton is " << clip_bottom << " clip_top is " << clip_top ); + + clip_row.AddNode( 0, SGGeod::fromDeg(-180.0, clip_bottom) ); + clip_row.AddNode( 0, SGGeod::fromDeg( 180.0, clip_bottom) ); + clip_row.AddNode( 0, SGGeod::fromDeg( 180.0, clip_top) ); + clip_row.AddNode( 0, SGGeod::fromDeg(-180.0, clip_top) ); + + clipped = tgPolygon::Intersect( clip_row, subject ); + if ( clipped.TotalNodes() > 0 ) { + ClipRow( clipped, b_clip.get_center_lat(), type ); + +#if 0 + { + char layer[32]; + char ds_name[64]; + sprintf(layer, "clipped_%d", strip_clip++ ); + sprintf(ds_name, "./stripped_%s", type.c_str() ); + + tgShapefile::FromPolygon( clipped, ds_name, layer, "poly" ); } - } - return; - } +#endif - // we have two or more rows left, split in half (along a - // horizontal dividing line) and recurse with each half - - // find mid point (integer math) - int mid = (dy + 1) / 2 - 1; - - // determine horizontal clip line - SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, mid); - double clip_line = b_clip.get_center_lat(); - if ( (clip_line >= -90.0 + SG_HALF_BUCKET_SPAN) - && (clip_line < 90.0 - SG_HALF_BUCKET_SPAN) ) - clip_line += SG_HALF_BUCKET_SPAN; - else if ( clip_line < -89.0 ) - clip_line = -89.0; - else if ( clip_line >= 89.0 ) - clip_line = 90.0; - else { - SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << clip_line ); - } - - { - // - // Crop bottom area (hopefully by putting this in it's own - // scope we can shorten the life of some really large data - // structures to reduce memory use) - // - - SG_LOG( SG_GENERAL, SG_DEBUG, "Generating bottom half (" << bb.getMin().getLatitudeDeg() << "-" << clip_line << ")" ); - - tgPolygon bottom, bottom_clip; - - bottom.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMin().getLatitudeDeg()) ); - bottom.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMin().getLatitudeDeg()) ); - bottom.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) ); - bottom.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) ); - - bottom_clip = tgPolygon::Intersect( subject, bottom ); - - if ( (bottom_clip.TotalNodes() > 0) && (bottom_clip.TotalNodes() != subject.TotalNodes()) ) { - Add( bottom_clip, type ); - } - } - - { - // - // Crop top area (hopefully by putting this in it's own scope - // we can shorten the life of some really large data - // structures to reduce memory use) - // - - SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half (" << clip_line << "-" << bb.getMax().getLatitudeDeg() << ")" ); - - tgPolygon top, top_clip; - - top.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) ); - top.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) ); - top.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMax().getLatitudeDeg()) ); - top.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMax().getLatitudeDeg()) ); - - top_clip = tgPolygon::Intersect( subject, top ); - - if ( (top_clip.TotalNodes() > 0) && (top_clip.TotalNodes() != subject.TotalNodes()) ) { - Add( top_clip, type ); } } } diff --git a/src/Lib/terragear/tg_chopper.hxx b/src/Lib/terragear/tg_chopper.hxx index af136366..3c526ff5 100644 --- a/src/Lib/terragear/tg_chopper.hxx +++ b/src/Lib/terragear/tg_chopper.hxx @@ -18,6 +18,7 @@ public: private: long int GenerateIndex( std::string path ); + void ClipRow( const tgPolygon& subject, const double& center_lat, const std::string& type ); tgPolygon Clip( const tgPolygon& subject, const std::string& type, SGBucket& b ); void Chop( const tgPolygon& subject, const std::string& type ); diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index 282a5287..35ba70c6 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -73,6 +73,23 @@ double tgContour::GetArea( void ) const return fabs(area * 0.5); } +bool tgContour::IsInside( const tgContour& inside, const tgContour& outside ) +{ + // first contour is inside second if the intersection of first with second is == first + // Intersection returns a polygon... + tgPolygon result; + bool isInside = false; + result = Intersect( inside, outside ); + + if ( result.Contours() == 1 ) { + if ( result.GetContour(0) == inside ) { + isInside = true; + } + } + + return isInside; +} + bool tgContour::RemoveCycles( const tgContour& subject, tgcontour_list& result ) { SG_LOG(SG_GENERAL, SG_DEBUG, "remove cycles : contour has " << subject.GetSize() << " points" ); @@ -116,6 +133,27 @@ bool tgContour::RemoveCycles( const tgContour& subject, tgcontour_list& result ) second.AddNode( subject.GetNode(n) ); } + // determine hole vs boundary + if ( IsInside( first, second ) ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "first contur is within second contour " ); + + // first contour is inside second : mark first contour as opposite of subject + first.SetHole( !subject.GetHole() ); + second.SetHole( subject.GetHole() ); + } else if ( IsInside( second, first ) ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "second contur is within first contour " ); + + // second contour is inside first : mark second contour as opposite of subject + first.SetHole( subject.GetHole() ); + second.SetHole( !subject.GetHole() ); + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, "conturs are (mostly) disjoint " ); + + // neither contour is inside - bots are the same as subject + first.SetHole( subject.GetHole() ); + second.SetHole( subject.GetHole() ); + } + SG_LOG(SG_GENERAL, SG_DEBUG, "remove first: size " << first.GetSize() ); first.SetHole( subject.GetHole() ); RemoveCycles( first, result ); @@ -460,6 +498,32 @@ tgPolygon tgContour::Diff( const tgContour& subject, tgPolygon& clip ) return result; } +tgPolygon tgContour::Intersect( const tgContour& subject, const tgContour& clip ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { + all_nodes.add( subject.GetNode(i) ); + } + + ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); + ClipperLib::Polygon clipper_clip = tgContour::ToClipper( clip ); + ClipperLib::Polygons clipper_result; + + ClipperLib::Clipper c; + c.Clear(); + c.AddPolygon(clipper_subject, ClipperLib::ptSubject); + c.AddPolygon(clipper_clip, ClipperLib::ptClip); + c.Execute(ClipperLib::ctIntersection, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + return result; +} + static bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, const std::vector& nodes, SGGeod& result, double bbEpsilon, double errEpsilon ) diff --git a/src/Lib/terragear/tg_contour.hxx b/src/Lib/terragear/tg_contour.hxx index 2063f02c..88a2bc90 100644 --- a/src/Lib/terragear/tg_contour.hxx +++ b/src/Lib/terragear/tg_contour.hxx @@ -77,6 +77,25 @@ public: double GetMinimumAngle( void ) const; double GetArea( void ) const; + bool operator==(const tgContour& other ) { + bool isEqual = true; + + if ( GetSize() == other.GetSize() ) + { + for (unsigned int i=0; i& nodes ); static bool FindColinearLine( const tgContour& subject, const SGGeod& node, SGGeod& start, SGGeod& end ); diff --git a/src/Prep/OGRDecode/ogr-decode.cxx b/src/Prep/OGRDecode/ogr-decode.cxx index aaef9238..eaa8cef3 100644 --- a/src/Prep/OGRDecode/ogr-decode.cxx +++ b/src/Prep/OGRDecode/ogr-decode.cxx @@ -63,6 +63,10 @@ string attribute_query; bool use_spatial_query=false; double spat_min_x, spat_min_y, spat_max_x, spat_max_y; int num_threads = 1; +bool save_shapefiles=false; +std::string ds_name="."; + +const double gSnap = 0.00000001; // approx 1 mm SGLockedQueue global_workQueue; @@ -197,6 +201,7 @@ void Decoder::processPolygon(OGRPolygon* poGeometry, const string& area_type ) // first add the outer ring tgPolygon shape = tgShapefile::ToPolygon( poGeometry ); + shape = tgPolygon::Simplify( shape ); if ( max_segment_length > 0 ) { shape = tgPolygon::SplitLongEdges( shape, max_segment_length ); @@ -608,6 +613,10 @@ int main( int argc, char **argv ) { num_threads=boost::thread::hardware_concurrency(); argv+=1; argc-=1; + } else if (!strcmp(argv[1],"--debug")) { + argv++; + argc--; + save_shapefiles=true; } else if (!strcmp(argv[1],"--help")) { usage(progname); } else {