- clipper fix for hole issue

- some missed code cleanup

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)

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.

src/Lib/Polygon/clipper.cpp

@@ -1,8 +1,8 @@
 /*******************************************************************************
 *                                                                              *
 * Author    :  Angus Johnson                                                   *
-* Version   :  4.8.9                                                           *
+* Version   :  4.9.1                                                           *
-* Date      :  25 September 2012                                               *
+* 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) ==
+    return Int128(e1.deltaY) * Int128(e2.deltaX) ==
-      Int128(e1.xtop - e1.xbot) * Int128(e2.ytop - e2.ybot);
+      Int128(e1.deltaX) * Int128(e2.deltaY);
-  else return (e1.ytop - e1.ybot)*(e2.xtop - e2.xbot) ==
+  else return e1.deltaY * e2.deltaX == e1.deltaX * e2.deltaY;
-      (e1.xtop - e1.xbot)*(e2.ytop - e2.ybot);
 }
 //------------------------------------------------------------------------------
 
@@ -539,8 +538,11 @@ double GetDx(const IntPoint pt1, const IntPoint pt2)
 
 void SetDx(TEdge &e)
 {
-  if (e.ybot == e.ytop) e.dx = HORIZONTAL;
+  e.deltaX = (e.xtop - e.xbot);
-  else e.dx = (double)(e.xtop - e.xbot) / (double)(e.ytop - e.ybot);
+  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) ?
+  return (val < 0) ? static_cast<long64>(val - 0.5) : static_cast<long64>(val + 0.5);
-    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;
+      if( !ePrev ) e = m_ActiveEdges;
-      else e = ePrior->nextInAEL;
+      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 &&
+      TEdge* ePrev = e->prevInAEL;
-        e->prevInAEL->xcurr == e->xbot && e->prevInAEL->ycurr == e->ybot &&
+      TEdge* eNext = e->nextInAEL;
-        SlopesEqual(IntPoint(e->xbot,e->ybot), IntPoint(e->xtop, e->ytop),
+      if (ePrev && ePrev->xcurr == e->xbot &&
-          IntPoint(e->xbot,e->ybot),
+        ePrev->ycurr == e->ybot && e->outIdx >= 0 &&
-          IntPoint(e->prevInAEL->xtop, e->prevInAEL->ytop), m_UseFullRange))
+        ePrev->outIdx >= 0 && ePrev->ycurr > ePrev->ytop &&
+        SlopesEqual(*e, *ePrev, m_UseFullRange))
       {
-        AddOutPt(e->prevInAEL, IntPoint(e->xbot, e->ybot));
+        AddOutPt(ePrev, IntPoint(e->xbot, e->ybot));
-        AddJoin(e, e->prevInAEL);
+        AddJoin(e, ePrev);
       }
-      else if (e->outIdx >= 0 && e->nextInAEL && e->nextInAEL->outIdx >= 0 &&
+      else if (eNext && eNext->xcurr == e->xbot &&
-        e->nextInAEL->ycurr > e->nextInAEL->ytop &&
+        eNext->ycurr == e->ybot && e->outIdx >= 0 &&
-        e->nextInAEL->ycurr <= e->nextInAEL->ybot &&
+        eNext->outIdx >= 0 && eNext->ycurr > eNext->ytop &&
-        e->nextInAEL->xcurr == e->xbot && e->nextInAEL->ycurr == e->ybot &&
+        SlopesEqual(*e, *eNext, m_UseFullRange))
-        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))
       {
-        AddOutPt(e->nextInAEL, IntPoint(e->xbot, e->ybot));
+        AddOutPt(eNext, IntPoint(e->xbot, e->ybot));
-        AddJoin(e, e->nextInAEL);
+        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); //nb: do this BEFORE testing orientation
-      FixupOutPolygon(*outRec1);
+      FixupOutPolygon(*outRec2); //    but AFTER calling PointIsVertex()
-      FixupOutPolygon(*outRec2);
 
+      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)

src/Lib/Polygon/clipper.hpp

@@ -1,8 +1,8 @@
 /*******************************************************************************
 *                                                                              *
 * Author    :  Angus Johnson                                                   *
-* Version   :  4.8.9                                                           *
+* Version   :  4.9.1                                                           *
-* Date      :  25 September 2012                                               *
+* 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);