diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index 65be2329..e875c79e 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -29,28 +29,9 @@ # include <config.h> #endif -//#include <simgear/compiler.h> - -//#include <iostream> -//#include <string> -//#include <vector> -//#include <algorithm> - -//#include <simgear/constants.h> -//#include <simgear/bucket/newbucket.hxx> #include <simgear/debug/logstream.hxx> -//#include <simgear/misc/sg_dir.hxx> -//#include <simgear/misc/sg_path.hxx> -//#include <simgear/math/sg_types.hxx> - -//#include <simgear/math/SGMath.hxx> -//#include <simgear/misc/sgstream.hxx> - - -//#include <boost/foreach.hpp> #include <Geometry/poly_support.hxx> -//#include <landcover/landcover.hxx> #include "tgconstruct.hxx" #include "usgs.hxx" @@ -58,76 +39,6 @@ using std::string; using std::vector; -//using namespace std; - -vector<string> load_dirs; -double nudge=0.0; - -#if 0 -// For each triangle assigned to the "default" area type, see if we -// can lookup a better land cover type from the 1km data structure. -static void fix_land_cover_assignments( TGConstruct& c ) { - SG_LOG(SG_GENERAL, SG_ALERT, "Fixing up default land cover types"); - // the list of node locations - TGTriNodes trinodes = c.get_tri_nodes(); - point_list geod_nodes = trinodes.get_node_list(); - - // the list of triangles (with area type attribute) - triele_list tri_elements = c.get_tri_elements(); - - // traverse the triangle element groups - SG_LOG(SG_GENERAL, SG_ALERT, " Total Nodes = " << geod_nodes.size()); - SG_LOG(SG_GENERAL, SG_ALERT, " Total triangles = " << tri_elements.size()); - for ( unsigned int i = 0; i < tri_elements.size(); ++i ) { - TGTriEle t = tri_elements[i]; - if ( t.get_attribute() == get_default_area_type() ) { - Point3D p1 = geod_nodes[t.get_n1()]; - Point3D p2 = geod_nodes[t.get_n2()]; - Point3D p3 = geod_nodes[t.get_n3()]; - - // offset by -quarter_cover_size because that is what we - // do for the coverage squares - AreaType a1 = get_area_type( c.get_cover(), - p1.x() - quarter_cover_size, - p1.y() - quarter_cover_size, - cover_size, cover_size ); - AreaType a2 = get_area_type( c.get_cover(), - p2.x() - quarter_cover_size, - p2.y() - quarter_cover_size, - cover_size, cover_size ); - AreaType a3 = get_area_type( c.get_cover(), - p3.x() - quarter_cover_size, - p3.y() - quarter_cover_size, - cover_size, cover_size ); - - // update the original triangle element attribute - AreaType new_area; - - // majority rules - if ( a1 == a2 ) { - new_area = a1; - } else if ( a1 == a3 ) { - new_area = a1; - } else if ( a2 == a3 ) { - new_area = a2; - } else { - // a different coverage for each vertex, just pick - // from the middle/average - Point3D average = ( p1 + p2 + p3 ) / 3.0; - //cout << " average triangle center = " << average; - new_area = get_area_type( c.get_cover(), - average.x() - quarter_cover_size, - average.y() - quarter_cover_size, - cover_size, cover_size ); - } - - //cout << " new attrib = " << get_area_name( new_area ) << endl; - c.set_tri_attribute( i, new_area ); - } - } -} -#endif - // display usage and exit static void usage( const string name ) { SG_LOG(SG_GENERAL, SG_ALERT, "Usage: " << name); @@ -161,12 +72,14 @@ int main(int argc, char **argv) { double ydist = -1; long tile_id = -1; + vector<string> load_dirs; + bool ignoreLandmass = false; + double nudge=0.0; + string debug_dir = "."; vector<string> debug_shape_defs; vector<string> debug_area_defs; - bool ignoreLandmass = false; - sglog().setLogLevels( SG_ALL, SG_INFO ); // Initialize shapefile support (for debugging) diff --git a/src/BuildTiles/Main/tgconstruct_landclass.cxx b/src/BuildTiles/Main/tgconstruct_landclass.cxx index 48973cf1..62021048 100644 --- a/src/BuildTiles/Main/tgconstruct_landclass.cxx +++ b/src/BuildTiles/Main/tgconstruct_landclass.cxx @@ -146,6 +146,71 @@ AreaType TGConstruct::get_landcover_type (const LandCover &cover, double xpos, d return get_default_area_type(); } +#if 0 +// For each triangle assigned to the "default" area type, see if we +// can lookup a better land cover type from the 1km data structure. +static void fix_land_cover_assignments( TGConstruct& c ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Fixing up default land cover types"); + // the list of node locations + TGTriNodes trinodes = c.get_tri_nodes(); + point_list geod_nodes = trinodes.get_node_list(); + + // the list of triangles (with area type attribute) + triele_list tri_elements = c.get_tri_elements(); + + // traverse the triangle element groups + SG_LOG(SG_GENERAL, SG_ALERT, " Total Nodes = " << geod_nodes.size()); + SG_LOG(SG_GENERAL, SG_ALERT, " Total triangles = " << tri_elements.size()); + for ( unsigned int i = 0; i < tri_elements.size(); ++i ) { + TGTriEle t = tri_elements[i]; + if ( t.get_attribute() == get_default_area_type() ) { + Point3D p1 = geod_nodes[t.get_n1()]; + Point3D p2 = geod_nodes[t.get_n2()]; + Point3D p3 = geod_nodes[t.get_n3()]; + + // offset by -quarter_cover_size because that is what we + // do for the coverage squares + AreaType a1 = get_area_type( c.get_cover(), + p1.x() - quarter_cover_size, + p1.y() - quarter_cover_size, + cover_size, cover_size ); + AreaType a2 = get_area_type( c.get_cover(), + p2.x() - quarter_cover_size, + p2.y() - quarter_cover_size, + cover_size, cover_size ); + AreaType a3 = get_area_type( c.get_cover(), + p3.x() - quarter_cover_size, + p3.y() - quarter_cover_size, + cover_size, cover_size ); + + // update the original triangle element attribute + AreaType new_area; + + // majority rules + if ( a1 == a2 ) { + new_area = a1; + } else if ( a1 == a3 ) { + new_area = a1; + } else if ( a2 == a3 ) { + new_area = a2; + } else { + // a different coverage for each vertex, just pick + // from the middle/average + Point3D average = ( p1 + p2 + p3 ) / 3.0; + //cout << " average triangle center = " << average; + new_area = get_area_type( c.get_cover(), + average.x() - quarter_cover_size, + average.y() - quarter_cover_size, + cover_size, cover_size ); + } + + //cout << " new attrib = " << get_area_name( new_area ) << endl; + c.set_tri_attribute( i, new_area ); + } + } +} +#endif + // Generate polygons from land-cover raster. Horizontally- or // vertically-adjacent polygons will be merged automatically. diff --git a/src/Lib/Polygon/clipper.cpp b/src/Lib/Polygon/clipper.cpp index 631973e4..15e8b783 100644 --- a/src/Lib/Polygon/clipper.cpp +++ b/src/Lib/Polygon/clipper.cpp @@ -1,8 +1,8 @@ /******************************************************************************* * * * Author : Angus Johnson * -* Version : 4.8.9 * -* Date : 25 September 2012 * +* Version : 4.9.1 * +* Date : 9 October 2012 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2012 * * * @@ -503,10 +503,9 @@ bool PointInPolygon(const IntPoint &pt, OutPt *pp, bool UseFullInt64Range) bool SlopesEqual(TEdge &e1, TEdge &e2, bool UseFullInt64Range) { if (UseFullInt64Range) - return Int128(e1.ytop - e1.ybot) * Int128(e2.xtop - e2.xbot) == - Int128(e1.xtop - e1.xbot) * Int128(e2.ytop - e2.ybot); - else return (e1.ytop - e1.ybot)*(e2.xtop - e2.xbot) == - (e1.xtop - e1.xbot)*(e2.ytop - e2.ybot); + return Int128(e1.deltaY) * Int128(e2.deltaX) == + Int128(e1.deltaX) * Int128(e2.deltaY); + else return e1.deltaY * e2.deltaX == e1.deltaX * e2.deltaY; } //------------------------------------------------------------------------------ @@ -539,8 +538,11 @@ double GetDx(const IntPoint pt1, const IntPoint pt2) void SetDx(TEdge &e) { - if (e.ybot == e.ytop) e.dx = HORIZONTAL; - else e.dx = (double)(e.xtop - e.xbot) / (double)(e.ytop - e.ybot); + e.deltaX = (e.xtop - e.xbot); + e.deltaY = (e.ytop - e.ybot); + + if (e.deltaY == 0) e.dx = HORIZONTAL; + else e.dx = (double)(e.deltaX) / (double)(e.deltaY); } //--------------------------------------------------------------------------- @@ -562,8 +564,7 @@ void SwapPolyIndexes(TEdge &edge1, TEdge &edge2) inline long64 Round(double val) { - return (val < 0) ? - static_cast<long64>(val - 0.5) : static_cast<long64>(val + 0.5); + return (val < 0) ? static_cast<long64>(val - 0.5) : static_cast<long64>(val + 0.5); } //------------------------------------------------------------------------------ @@ -2135,12 +2136,13 @@ void Clipper::AddOutPt(TEdge *e, const IntPoint &pt) { //check for 'rounding' artefacts ... if (outRec->sides == esNeither && pt.Y == op->pt.Y) + { if (ToFront) { if (pt.X == op->pt.X +1) return; //ie wrong side of bottomPt } else if (pt.X == op->pt.X -1) return; //ie wrong side of bottomPt - + } outRec->sides = (EdgeSide)(outRec->sides | e->side); if (outRec->sides == esBoth) { @@ -2641,10 +2643,10 @@ void Clipper::ProcessEdgesAtTopOfScanbeam(const long64 topY) if( IsMaxima(e, topY) && !NEAR_EQUAL(GetMaximaPair(e)->dx, HORIZONTAL) ) { //'e' might be removed from AEL, as may any following edges so ... - TEdge* ePrior = e->prevInAEL; + TEdge* ePrev = e->prevInAEL; DoMaxima(e, topY); - if( !ePrior ) e = m_ActiveEdges; - else e = ePrior->nextInAEL; + if( !ePrev ) e = m_ActiveEdges; + else e = ePrev->nextInAEL; } else { @@ -2693,25 +2695,23 @@ void Clipper::ProcessEdgesAtTopOfScanbeam(const long64 topY) UpdateEdgeIntoAEL(e); //if output polygons share an edge, they'll need joining later ... - if (e->outIdx >= 0 && e->prevInAEL && e->prevInAEL->outIdx >= 0 && - e->prevInAEL->xcurr == e->xbot && e->prevInAEL->ycurr == e->ybot && - SlopesEqual(IntPoint(e->xbot,e->ybot), IntPoint(e->xtop, e->ytop), - IntPoint(e->xbot,e->ybot), - IntPoint(e->prevInAEL->xtop, e->prevInAEL->ytop), m_UseFullRange)) + TEdge* ePrev = e->prevInAEL; + TEdge* eNext = e->nextInAEL; + if (ePrev && ePrev->xcurr == e->xbot && + ePrev->ycurr == e->ybot && e->outIdx >= 0 && + ePrev->outIdx >= 0 && ePrev->ycurr > ePrev->ytop && + SlopesEqual(*e, *ePrev, m_UseFullRange)) { - AddOutPt(e->prevInAEL, IntPoint(e->xbot, e->ybot)); - AddJoin(e, e->prevInAEL); + AddOutPt(ePrev, IntPoint(e->xbot, e->ybot)); + AddJoin(e, ePrev); } - else if (e->outIdx >= 0 && e->nextInAEL && e->nextInAEL->outIdx >= 0 && - e->nextInAEL->ycurr > e->nextInAEL->ytop && - e->nextInAEL->ycurr <= e->nextInAEL->ybot && - e->nextInAEL->xcurr == e->xbot && e->nextInAEL->ycurr == e->ybot && - SlopesEqual(IntPoint(e->xbot,e->ybot), IntPoint(e->xtop, e->ytop), - IntPoint(e->xbot,e->ybot), - IntPoint(e->nextInAEL->xtop, e->nextInAEL->ytop), m_UseFullRange)) + else if (eNext && eNext->xcurr == e->xbot && + eNext->ycurr == e->ybot && e->outIdx >= 0 && + eNext->outIdx >= 0 && eNext->ycurr > eNext->ytop && + SlopesEqual(*e, *eNext, m_UseFullRange)) { - AddOutPt(e->nextInAEL, IntPoint(e->xbot, e->ybot)); - AddJoin(e, e->nextInAEL); + AddOutPt(eNext, IntPoint(e->xbot, e->ybot)); + AddJoin(e, eNext); } } e = e->nextInAEL; @@ -2940,30 +2940,6 @@ void Clipper::DoBothEdges(TEdge *edge1, TEdge *edge2, const IntPoint &pt) } //---------------------------------------------------------------------- -void Clipper::CheckHoleLinkages1(OutRec *outRec1, OutRec *outRec2) -{ - //when a polygon is split into 2 polygons, make sure any holes the original - //polygon contained link to the correct polygon ... - for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); ++i) - { - OutRec *orec = m_PolyOuts[i]; - if (orec->isHole && orec->bottomPt && orec->FirstLeft == outRec1 && - !PointInPolygon(orec->bottomPt->pt, outRec1->pts, m_UseFullRange)) - orec->FirstLeft = outRec2; - } -} -//---------------------------------------------------------------------- - -void Clipper::CheckHoleLinkages2(OutRec *outRec1, OutRec *outRec2) -{ - //if a hole is owned by outRec2 then make it owned by outRec1 ... - for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); ++i) - if (m_PolyOuts[i]->isHole && m_PolyOuts[i]->bottomPt && - m_PolyOuts[i]->FirstLeft == outRec2) - m_PolyOuts[i]->FirstLeft = outRec1; -} -//---------------------------------------------------------------------- - void Clipper::JoinCommonEdges(bool fixHoleLinkages) { for (JoinList::size_type i = 0; i < m_Joins.size(); i++) @@ -3049,32 +3025,25 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages) outRec2->bottomPt = outRec2->pts; outRec2->bottomPt->idx = outRec2->idx; + int K = 0; if (PointInPolygon(outRec2->pts->pt, outRec1->pts, m_UseFullRange)) { //outRec2 is contained by outRec1 ... + K = 1; outRec2->isHole = !outRec1->isHole; outRec2->FirstLeft = outRec1; - if (outRec2->isHole == - (m_ReverseOutput ^ Orientation(outRec2, m_UseFullRange))) - ReversePolyPtLinks(*outRec2->pts); } else if (PointInPolygon(outRec1->pts->pt, outRec2->pts, m_UseFullRange)) { //outRec1 is contained by outRec2 ... + K = 2; outRec2->isHole = outRec1->isHole; outRec1->isHole = !outRec2->isHole; outRec2->FirstLeft = outRec1->FirstLeft; outRec1->FirstLeft = outRec2; - if (outRec1->isHole == - (m_ReverseOutput ^ Orientation(outRec1, m_UseFullRange))) - ReversePolyPtLinks(*outRec1->pts); - //make sure any contained holes now link to the correct polygon ... - if (fixHoleLinkages) CheckHoleLinkages1(outRec1, outRec2); } else { outRec2->isHole = outRec1->isHole; outRec2->FirstLeft = outRec1->FirstLeft; - //make sure any contained holes now link to the correct polygon ... - if (fixHoleLinkages) CheckHoleLinkages1(outRec1, outRec2); } //now fixup any subsequent joins that match this polygon @@ -3087,10 +3056,47 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages) j2->poly2Idx = j->poly2Idx; } - //now cleanup redundant edges too ... - FixupOutPolygon(*outRec1); - FixupOutPolygon(*outRec2); + FixupOutPolygon(*outRec1); //nb: do this BEFORE testing orientation + FixupOutPolygon(*outRec2); // but AFTER calling PointIsVertex() + switch( K ) { + case 1: + { + if (outRec2->isHole == + (m_ReverseOutput ^ Orientation(outRec2, m_UseFullRange))) + ReversePolyPtLinks(*outRec2->pts); + break; + } + case 2: + { + if (outRec1->isHole == + (m_ReverseOutput ^ Orientation(outRec1, m_UseFullRange))) + ReversePolyPtLinks(*outRec1->pts); + //make sure any contained holes now link to the correct polygon ... + if (fixHoleLinkages && outRec1->isHole) + for (PolyOutList::size_type k = 0; k < m_PolyOuts.size(); ++k) + { + OutRec *orec = m_PolyOuts[k]; + if (orec->isHole && orec->bottomPt && orec->FirstLeft == outRec1) + orec->FirstLeft = outRec2; + } + break; + } + default: + { + //make sure any contained holes now link to the correct polygon ... + if (fixHoleLinkages) + for (PolyOutList::size_type k = 0; k < m_PolyOuts.size(); ++k) + { + OutRec *orec = m_PolyOuts[k]; + if (orec->isHole && orec->bottomPt && orec->FirstLeft == outRec1 && + !PointInPolygon(orec->bottomPt->pt, outRec1->pts, m_UseFullRange)) + orec->FirstLeft = outRec2; + } + + } + } + if (Orientation(outRec1, m_UseFullRange) != (Area(*outRec1, m_UseFullRange) > 0)) DisposeBottomPt(*outRec1); if (Orientation(outRec2, m_UseFullRange) != (Area(*outRec2, m_UseFullRange) > 0)) @@ -3101,9 +3107,14 @@ void Clipper::JoinCommonEdges(bool fixHoleLinkages) //joined 2 polygons together ... //make sure any holes contained by outRec2 now link to outRec1 ... - if (fixHoleLinkages) CheckHoleLinkages2(outRec1, outRec2); + if (fixHoleLinkages) + for (PolyOutList::size_type k = 0; k < m_PolyOuts.size(); ++k) + if (m_PolyOuts[k]->isHole && m_PolyOuts[k]->bottomPt && + m_PolyOuts[k]->FirstLeft == outRec2) + m_PolyOuts[k]->FirstLeft = outRec1; - //now cleanup redundant edges too ... + + //and cleanup redundant edges too ... FixupOutPolygon(*outRec1); if (outRec1->pts) diff --git a/src/Lib/Polygon/clipper.hpp b/src/Lib/Polygon/clipper.hpp index 351ac79c..a3427f74 100644 --- a/src/Lib/Polygon/clipper.hpp +++ b/src/Lib/Polygon/clipper.hpp @@ -1,8 +1,8 @@ /******************************************************************************* * * * Author : Angus Johnson * -* Version : 4.8.9 * -* Date : 25 September 2012 * +* Version : 4.9.1 * +* Date : 9 October 2012 * * Website : http://www.angusj.com * * Copyright : Angus Johnson 2010-2012 * * * @@ -98,6 +98,8 @@ struct TEdge { long64 xtop; long64 ytop; double dx; + long64 deltaX; + long64 deltaY; long64 tmpX; PolyType polyType; EdgeSide side; @@ -276,8 +278,6 @@ private: void FixupOutPolygon(OutRec &outRec); bool IsHole(TEdge *e); void FixHoleLinkage(OutRec *outRec); - void CheckHoleLinkages1(OutRec *outRec1, OutRec *outRec2); - void CheckHoleLinkages2(OutRec *outRec1, OutRec *outRec2); void AddJoin(TEdge *e1, TEdge *e2, int e1OutIdx = -1, int e2OutIdx = -1); void ClearJoins(); void AddHorzJoin(TEdge *e, int idx);